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ABSTRACT 
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The  unsteady,  compressible,  Reynolds-averaged  Navier-Stokes 
equations  were  solved  for  the  flow  field  about  an  external  compression 
axi-symmetric  inlet  with  a  length  to  diameter  ratio,  L/D  *  15.88,  at 
Mach  2.0  and  a  Reynolds  number  based  on  diameter,  Re^  ■  2.36x10**,  operat¬ 
ing  in  the  near-critical  and  subcritical  flow  regimes.  The  near-critical 
solution  reached  a  stable  steady  state  while  the  subcritical  solutions 
attained  an  unstable  bounded  oscillatory  state,  characterized  by  large 
amplitude  pressure  oscillations  and  traveling  shock  waves.  This  phenomenon 
is  the  result  of  a  shear  layer  instability  combined  with  a  closed-loop 
feedback  of  reflected  disturbances  and  the  naturally  occurring  self- 
sustained  oscillations  are  commonly  known  as  buzz.  Numerical  results 
are  given  in  terms  of  Mach  contours,  velocity  plots,  pressure-time  traces 
at  selected  stations,  as  well  as  mass  flux  and  other  mass-averaged 
quantities  along  the  duct  length.  Comparison  with  experiment  is  also 
given. 
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Section  I 


INTRODUCTION 

The  fluid  dynamic  behavior  of  a  supersonic  inlet,  operating  in  the 
critical  and  subcrltlcal  flow  regimes  has  been  of  interest  since  the 
earliest  experiments  with  supersonic  diffusers.  These  flow  conditions 
correspond  to  a  normal  shock  location  at  or  before  the  diffuser  throat 
(minimum  cross  sectional  area)  near  the  inlet  entrance.  The  critical 
flow  regime  is  the  point  of  optimum  diffuser  efficiency  and  results 
in  maximum  mass  flow,  pressure  recovery,  and  stream  thrust.  Unfortunately, 
the  critical  and,  more  typically,  the  subcritical  flow  regimes  are 
characterized  by  the  onset  of  violent  flow  oscillations  within  the 
inlet  and  rapid  movement  of  the  normal  or  bow  shock  on  the  centerbody. 

This  motion  is  termed  buzz  and  is  frequently  violent  enough  to  result  in 
structural  damage,  engine  surge  (jet  engine),  combustion  flaaeout,  or 
non-recoverable  thrust  loss  (ramjet).  Although  the  buzz  phenomenon  has 
been  studied  for  almost  forty  years  it  remains  poorly  understood.  To 
this  day,  there  is  no  analytic  method  for  predicting  the  amplitude  or 
frequency  of  the  oscillations  or  even  a  certain  method  for  predicting  buzz 
onset. 

With  the  advent  of  "super  computers"  such  as  the  CRAY  1-S  and  the 
rapid  maturity  of  computational  fluid  dynamics  as  a  discipline,  it  has 
recently  become  possible  to  calculate  numerical  solutions  in  these  flow 
regimes  from  the  fundamental  equations  of  fluid  mechanics  -  the  Navier- 
Stokes  equations.  Although  this  approach  is  certainly  not  without  its 
own  difficulties,  it  does  offer  the  possibility  of  complete  information 
without  regard  to  instrumentation  and  the  potential  for  new  insight 
into  the  problem.  The  purpose  of  this  report  is  to  document  a  numerical 


study  of  the  near-critical  and  subcritical  regimes  of  an  external 
compression  axi-symmetric  inlet  tested  experimentally  at  the  University 
of  Tokyo  by  Nagashlma,  Obokata,  and  Asanuma  (Ref  1)  in  1972. 


1.1  Operating  Characteristics  of  Supersonic  Inlets 


The  operating  characteristics  of  a  supersonic  inlet  can  be  classi¬ 
fied  as  either  supercritical,  critical  or  subcritical  (Refs  2,3).  The 
three  regimes  are  illustrated  in  Figure  (1)  for  an  external  compression 
inlet  at  the  design  Mach  number.  In  this  type  of  inlet,  the  minimum 
entrance  area,  Ac>  occurs  at  the  cowl  lip.  The  design  Mach  number 
corresponds  to  a  flight  condition  in  which  the  conical  shock  impinges 
on  the  tip  of  the  cowl  lip.  For  a  given  flight  condition  (Mach  number 
and  angle  of  attack) ,  the  operating  condition  of  the  diffuser  is  de¬ 
termined  by  the  diffuser  back  pressure  required  to  match  that  dictated 
by  the  downstream  component.  This  downstream  component  could  be  a 
throat,  a  combustor,  or  a  compressor  face.  For  the  case  at  hand,  the 
downstream  component  is  a  throat  and  the  operating  condition  is  deter¬ 
mined  by  the  requirement  that  the  mass  flow  passed  by  the  downstream 
throat  must  equal  that  captured  by  the  inlet 
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and  all  properties  are  mass-averaged  over  the  station  cross  sectional 
area.  This  leads  to  the  relation 
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expelled  shock  wave 


Fig. 


Operating  Regimes  of  a  Supersonic  External 
Compression  Inlet  (from  Ref  3) 


In  the  supercritical  regime,  the  downstream  throat  may  be  assumed  to  be 
choked  for  the  throat  areas  considered.  The  mass  flow  captured  by  the 
inlet  is  a  maximum  and  is  determined  by  the  geometry  and  flight  condi¬ 
tion  (A  ,pt  ,  X  specified).  As  the  throat  area,  A*,  is  changed,  the 
c 

downstream  total  pressure,  p  ,  is  then  determined  by  Equation  (1.2). 

This  serves  to  position  the  shock  system  within  the  diffuser  since  the 
majority  of  the  total  pressure  loss  is  shock  induced.  As  the  downstream 
throat  is  further  closed,  a  higher  p  is  required.  This  can  only  be 
produced  if  the  shock  is  forced  upstream  in  the  diffuser  where  it  occurs 
at  a  lower  averaged  Mach  number  and  a  smaller  total  pressure  loss.  If  the 
throat  area  is  sufficiently  reduced  the  normal  shock  will  be  positioned 
at  the  cowl  lip.  As  mentioned  previously,  this  critical  regime  is  the 
point  where  maximum  mass  flow  and  highest  pressure  recovery  are  coin¬ 
cident.  Any  further  decrease  in  throat  area  will  force  the  normal  shock 
out  of  the  diffuser  to  a  position  on  the  centerbody  so  that  mass  flow 
may  be  spilled  to  a  level  which  can  be  passed  by  the  downstream  throat. 
This  is  the  subcrltlcal  regime  and  is  where  buzz  typically  occurs. 

In  the  past,  inlets  have  been  designed  to  avoid  the  unstable 
flow  regime  through  the  use  of  such  devices  as  by-pass  doors,  bound¬ 
ary  layer  bleed,  or  fuel  flow  limiters.  However,  the  requirement 
for  a  stable  operating  regime  often  imposes  a  substantial  penalty 
on  the  inlet  in  terms  of  diffuser  efficiency  (total  pressure  recovery) 
or  additive  drag.  This  impacts  total  flight  vehicle  performance  but 
most  significantly  it  may  impose  bounds  in  the  flight  envelope  for  off- 
design  and  maneuvering  conditions.  A  more  thorough  understanding  of 


buzz  is  a  necessary  first  step  toward  the  design  of  future  inlets 
less  susceptible  to  its  adverse  effects. 


1.2  Review  of  Previous  Investigations  into  Inlet  Buzz 

Oswatitsch  (Ref  4)  first  observed  inlet  buzz  during  the  course  of 
experimental  work  in  1944  but  did  not  explain  its  occurrence. 

Ferri  and  Nucci  (Ref  5)  tested  a  variety  of  external  compression 
diffuser  geometries.  They  identified  the  cause  of  the  instability  as 
the  entrance  of  a  vortex  sheet  (free  shear  layer)  into  the  cowl.  The 
vortex  sheet  resulted  from  (1)  the  Interaction  of  the  oblique  conical 
and  the  strong  bow  shock  or  (2)  a  lambda  shock  at  the  base  of  the 
centerbody  resulting  from  the  shock-induced  boundary  layer  separation. 
In  their  view,  the  vortex  sheet  tended  to  cause  the  boundary  layer 
on  the  inner  cowl  wall  to  separate  which  then  choked  the  flow.  Pres¬ 
sure  waves  propagated  forward  and  drove  the  shock  further  outward.  The 
subsequent  reduction  in  mass  flow  decreased  the  back  pressure,  the  bow 
shock  retreated  and  the  vortex  sheet  was  expelled.  Once  initiated, 
the  process  continued  in  a  periodic  fashion.  Their  theory  is  incom¬ 
plete  since  stable  flow  has  been  observed  with  vortex  sheet  impinge¬ 
ment  on  the  cowl. 

Sterbentz  and  Eward  (Ref  6)  and  Sterbentz  and  Davids  (Ref  7)  used 
a  Helmholz  resonator  analogy  in  which  the  inlet  mass  flow  was  modeled 
as  a  slug  of  air  bounded  downstream  by  a  weightless  spring  to  predict 
resonant  frequency  and  amplitude.  This  led  to  an  instability  require¬ 
ment  that  the  diffuser  pressure  recovery  curve  as  a  function  of  mass 
flow  should  have  a  positive  slope  in  the  subcritlcal  region.  This  hy¬ 


pothesis  has  been  shown  to  be  false.  Further,  the  theory  cannot  predict 


the  onset  of  buzz  and  does  not  model  the  actual  physical  processes. 

Trimpi  (Ref  8)  advanced  a  modified  one-dimensional  characteristic 
theory  in  which  the  buzz  cycle  was  traced  through  the  inlet  as  a  series 
of  compression  and  expansion  waves  propagating  and  reflecting  from 
the  bow  shock  and  the  choked  nozzle  throat.  His  theory  did  not  predict 
the  onset  of  buzz  and  in  fact  required  the  use  of  experimental  data 
to  start  the  calculation  when  the  shock  was  in  motion  ahead  of  the  cowl. 
The  Helmholz  resonator  and  wave  propagation  theories  are  not  compatible. 
Trimpi,  in  a  later  experimental  study  (Ref  9),  proposed  a  prediction 
method  for  buzz  amplitude  and  further  suggested  a  qualitative  method 
in  which  stability  was  related  to  the  instantaneous  ratio  of  entropy 
increase  to  mass  flow  decrease.  He  observed  the  occurrence  of  low 
frequency  buzz,  which  agreed  well  with  his  theory,  and  of  an  intermit¬ 
tent  high  frequency  buzz  whose  origin  he  postulated  to  be  localized 
near  the  inlet  entrance. 

Dailey  (Ref  10)  also  observed  experimentally  both  high  and  low 
frequency  oscillations.  For  high  frequency  oscillation  which  occurred 
at  low  mass  flows,  he  suggested  vortex  shedding  at  the  cowl  lip.  For 
low  frequency  buzz,  characteristic  at  higher  subcrltical  mass  flows, 

Dailey  postulated  that  the  instability  was  initiated  within  the  sep¬ 
arated  boundary  layer  on  the  centerbody.  This  is  in  contrast  with  the 
theory  of  Ferri  in  which  instability  resulted  from  vortex  sheet  im¬ 
pingement  on  the  cowl.  According  to  Dailey,  the  inlet  flow  was  completely 
blocked  due  to  massive  separation  within  the  inlet.  Upstream  of  the 
blockage,  mass  accumulated  and  forced  the  bow  shock  upstream  until  flow 
was  spilled  from  the  cowl  lip.  Downstream  of  the  blockage,  back  pres¬ 
sure  decreased  and  the  flow  collapsed  to  restart  the  cycle.  In  Dailey's 


view,  the  buzz  cycle  was  governed  by  a  relaxation  process  determined 
by  the  time  required  for  mass  to  fill  up  and  spill  from  the  diffuser. 
Trlmpi  had  previously  described  this  process  in  terms  of  wave  propaga¬ 
tion  and  Refs  (11)  and  (12)  indicate  the  sharp  disagreement  between 
the  theories  of  Trimpi  and  Dailey. 

In  a  subsequent  study.  Brown  (Ref  13)  tested  a  series  of  inlet 
types.  Results  were  similar  to  those  of  Dailey  (Ref  10).  Brown 
advanced  the  hypothesis  of  an  edge  tone  phenomenon  to  explan  the  high 
frequency  oscillation  and  the  sudden  shifts  in  frequency. 

Mlrels  (Ref  14)  performed  a  one  dimensional  analysis  of  small  per¬ 
turbations  in  which  the  mechanism  leading  to  buzz  is  the  amplifica¬ 
tion  of  acoustic  waves  due  to  successive  reflections  in  the  combus¬ 
tion  or  plenum  chamber.  He  concluded  the  flow  would  be  unstable  when 
the  real  part  of  the  acoustic  impedance  of  the  inlet  is  greater  than 
a  term  of  the  same  order  of  magnitude  as  the  combustion  chamber  Mach 
number . 

Previously,  Stoolman  (Ref  15)  had  considered  the  stability  of  a 
normal  shock  diffuser.  The  mechanism  for  stability  considered,  was 
again  the  amplification  ratio  of  reflected  to  incident  waves.  He  con¬ 
cluded  that  the  diffuser  was  unstable  in  the  absence  of  viscous  dis¬ 
sipation.  His  conclusion  was  contradicted  by  Chang  and  Hsu  (Ref  16) 
who  reported  that  acoustic  Impedance  was  always  negative  and  the  flow 
always  stable  in  the  absence  of  viscous  dissipation.  In  a  later 
paper  (Ref  17) ,  solutions  to  one-dimensional  quasi-steady,  linearized 
equations  for  viscous  compressible  flow  indicated  that  Instability  was 
possible,  due  to  wave  amplification  at  the  inlet  entrance  when  entropy 
production  due  to  viscous  dissipation  was  considered. 


Kowalewicz  (Ref  18)  considered  the  stability  of  a  normal  shock 
diffuser  with  the  extension  of  accounting  for  mass  flow  spillage  and 
divergence  in  diffuser  cross-sectional  area.  The  approach  was  to 
linearize  the  governing  equations,  assume  harmonic  solutions,  and  solve 
the  eigenvalue  problem  to  predict  natural  frequency  and  damping.  Since 
viscosity  was  not  considered,  for  high  subcritical  mass  flows  (low 
spillage,  where  the  theory  was  valid),  the  flow  was  found  to  be  stable, 
in  agreement  with  the  results  of  Chang  and  Hsu. 

In  experimental  testing  of  2-D  inlets  of  the  Concorde  type,  Fisher, 
Neale,  and  Brooks  (Ref  19)  identified  two  distinct  amplitudes  of  buzz 
of  roughly  the  same  frequency  -  little  buzz  and  big  buzz.  In  both 
cases  the  triggering  mechanism  was  the  entrance  of  a  shear  zone  with 
a  total  pressure  gradient  into  the  diffuser  (Ferri  condition) .  In¬ 
stability  was  dependent  upon  the  magnitude  of  the  total  pressure  change 
and  its  gradient,  as  well  as  its  proximity  to  the  cowl  lip.  Specifi¬ 
cally,  little  buzz  was  associated  with  the  intersection  of  the  normal 
shock  with  the  ramp  hinge  shock  and  isentropic  fan.  Big  buzz  was 
found  to  be  initiated  by  the  intersection  of  the  normal  shock  with 
the  initial  wedge  shock.  The  strength  of  the  shear  zone  total  pressure 
gradient  was  then  sufficient  to  cause  downstream  boundary  layer  sep¬ 
aration  which  resulted  in  inlet  throat  choking  and  large  amplitude 
oscillation  in  the  manner  described  by  Dailey. 

Steward  and  Fisher  (Ref  20)  studied  the  effects  of  freestream 
turbulence,  Reynolds  number,  and  surface  roughness,  on  axi-symmetric 
inlets  at  Mach  numbers  near  two.  They  observed  the  Ferri  instability 
at  centerbody  positions  in  which  the  shock  impinged  on  the  cowl  lip. 

For  other  positions,  they  found  the  instability  resulted  from  centerbody 


boundary  layer  instability.  At  higher  subcritlcal  mass  flows  this 
type  of  instability  resulted  in  a  very  mild  and  erratic  oscillation 
termed  flutter.  The  much  larger  amplitude  oscillation,  buss,  was 
encountered  if  a  turbulent  boundary  layer  trip  was  placed  at  the  center- 
body  tip  or  if  the  mass  flow  were  further  reduced. 

Nagashima,  Obokata  and  Ananuma  (Ref  1)  experimentally  studied  in¬ 
stability  in  axlsymmetrlc  external  compression  inlets.  They  found 
that  both  high  and  low  frequency  buzz,  as  measured  by  spectral  analyser, 
have  definite  spikes  as  predicted  by  an  acoustic  model.  This  report 
was  chosen  as  the  basis  for  numerical  comparison  and  is  described  in 
detail  in  Section  IV. 

Sajben,  et.al.  (Refs  21,  22,  23,  24)  appear  to  have  conducted  the 

most  comprehensive  experimental  investigation  of  diffuser  instability. 

Two  different  nominally  two  dimensional  supercritically  operating 

diffusers  with  area  ratios  of  A  /A  -  1.52  and  2.37  were  considered.  A 

e  c 

variety  of  pressure  ratios  were  tested,  corresponding  to  different  up¬ 
stream  normal  shock  Mach  numbers.  The  smaller  area  ratio  diffuser  was 
found  to  be  relatively  quiet  with  naturally  occurring  pressure  dis¬ 
turbances  on  the  order  of  AP/Pfc  *  .06,  while  the  larger  area  ratio 
model  displayed  disturbances  of  AP/Pt  -  .11  -*■  .15  and  shock  wave 
movement  equal  to  one  throat  height.  The  characteristic  length  scale 
for  the  weaker  shock  conditions  was  the  duct  length  and  the  disturbance 
propagation  was  acoustic  in  origin.  For  the  stronger  shock  conditions, 
the  characteristic  length  scale  corresponded  with  the  inviscid  core 
length  and  the  disturbance  propagation  was  not  acoustic.  Forced 
disturbances  at  the  diffuser  exit  were  also  investigated  and  resonance 


was  not  found 


1.3  Computational  Fluid  Dynamics  Applied  to  Internal  Flows  and  Flows 
With  Self-Sustained  Oscillations 
Internal  Flows 

The  field  of  computational  aerodynamics  has  continued  its  rapid 
evolutionary  growth  as  a  major  tool  in  aerodynamics  research.  Early 


computational  approaches  to  internal  flow  problems  consisted  of  the 
method  of  characteristics  or  potential  flow  methods  for  the  inviscid 
flow  coupled  with  a  boundary  layer  method  for  the  viscous  flow  field. 
However,  the  unsteady  Reynolds-averaged  Navier-Stokes  equations  must 
be  solved  for  many  "real  life"  diffusers  because  of  the  presence  of 
strong  shock-boundary  layer  interactions,  subsonic  flow,  and  boundary 
layer  separation  due  to  the  adverse  pressure  gradients  inherent  in 
diffusers.  Restricting  attention  to  the  Navier-Stokes  equations, 

Knight  (Ref  25)  gave  the  first  solution  to  a  realistic  high  Reynolds 
number  inlet.  Numerical  solutions  were  obtained  for  the  two-dimen¬ 
sional  steady  flow  about  two  hypersonic  inlet  geometries.  More  recently, 
Knight  (Refs  26,  27)  considered  the  supersonic  flow  associated  with  a 
Mach  2.5  inlet  with  the  extension  of  accounting  for  boundary  layer 
bleed  and  improving  algorithm  efficiency.  Currently,  work  on  three 
dimensional  Navier-Stokes  calculations  of  inlet  flows  is  still  in  pro¬ 
gress  (Ref  28) .  As  an  alternative  to  the  considerable  cost  of  solutions 
to  the  full  Navier-Stokes  equations,  parabolized  procedures  have  also 
been  developed,  in  which  streamwise  diffusive  terms  are  dropped,  allowing 
spatial  marching  in  the  primary  flow  direction.  Reference  (29)  describes 
a  procedure  which  has  been  applied  to  supersonic  diffuser  flows  with 
adverse  pressure  gradients  so  long  as  streamwise  separation  was  not 
encountered.  References  (30,  31)  summarize  the  status  of  current  compute- 


tiona1  methods  applied  to  internal  aerodynamics. 


Flows  Exhibiting  Self-Sustained  Oscillations 


Steady  state  solutions  to  the  Navier-Stokes  equations  are  nor¬ 
mally  obtained  from  the  unsteady  equations  by  a  relaxation  process  to 
an  asymptotic  steady  state  condition.  It  was  realized,  however,  that 
these  equations  should  also  provide  accurate  temporal  resolution  of 
unsteady  flows  as  well.  An  important  subset  of  unsteady  flows  includes 
those,  with  no  external  forcing  function,  for  which  a  steady  state 
solution  does  not  exist.  These  flows  exhibit  bounded  self-sustained 
oscillations  which  are  periodic  or  quasi-periodic.  Levy  (Ref  32)  was  the 
first  to  demonstrate  that  current  Navier-Stokes  algorithms  were  capable 
of  reproducing  the  time  dependent  aspects  of  the  mean  flow  in  an 
unsteady  turbulent  flow  field.  Levy  considered  the  flow  about  a  thick 
airfoil  in  the  transonic  regime.  The  flow  was  found  to  be  steady 
except  in  a  critical  Mach  number  regime  where  the  shock  on  the  upper 
and  lower  surfaces  displayed  an  asymmetric  oscillation.  The  flow 
oscillated  between  the  trailing  edge  separation  observed  at  lower  Mach 
numbers  and  a  shock  induced  separation  at  higher  Mach  numbers.  Good 
agreement  with  experiment  was  found  but  inadequacies  in  the  turbulence 
model  for  the  separated  region  were  noted. 

Steger  and  Bailey  (Ref  33)  successfully  simulated  the  fluid  dy¬ 
namically  driven  flutter  or  buzz  of  the  aileron  of  a  P-80  airfoil  at  the 
flight  condition  where  it  had  been  observed  experimentally  some  30  years 
earlier.  The  flow  was  steady  at  lower  Mach  numbers  and  unsteady  at  higher 
Mach  numbers,  in  accord  with  experiment.  Inviscid  calculations  were 


found  to  decay  or  diverge,  depending  upon  the  Mach  number,  indicating 


the  crucial  role  played  by  viscosity  in  establishing  a  repeatable,  bounded 
limit  cycle. 

Shang  and  Hankey  have  considered  several  flows  which  exhibit  self- 
sustained  oscillations.  The  pressure  oscillations  produced  by  the 
shear  layer  instability  in  the  supersonic  flow  over  an  open  cavity  were 
computed  in  Ref  (34) .  Linear  stability  theory  was  also  used  to  explain 
the  resonant  phenomenon.  In  Ref  (35) ,  numerical  solutions  were  obtained 
for  two  different  spike  lengths  for  the  axi-symmetric  flow  over  spike- 
tipped  bodies  in  a  Mach  three  laminar  flow  regime.  The  shorter  length 
body  resulted  in  a  stable  flow  in  accord  with  the  experiment  and  linear 
stability  theory.  The  longer  length  body  produced  a  large  amplitude 
high  frequency  oscillation  in  excellent  agreement  with  the  experiment 
in  both  frequency  and  amplitude  of  the  observed  oscillations.  Shang 
(Ref  36)  also  considered  the  oscillatory  flow  about  a  cylinder  involv¬ 
ing  the  periodic  shedding  of  vorticity.  The  numerical  simulation  t^ain 
was  successful  in  reproducing  the  dominant  fluid  dynamic  mechanisms. 

Very  recently,  Liou  and  Coakley  (Refs  37,  38)  computed  several 
of  the  cases  investigated  experimentally  by  Sajben  for  both  forced  and 
self-sustained  oscillations.  The  forced  oscillations  computed  were 
predicted  well  in  terms  of  amplitude  and  phase  angle  although  shock 
displacement  was  overpredicted  in  some  cases.  For  the  self-sustained 
oscillations,  good  agreement  in  predicted  frequency  was  obtained  at 
moderate  and  weak  shock  pressure  ratios.  At  the  strong  shock  pressure 
ratio,  the  computed  frequency  was  1.75  times  the  experimental  result. 

The  authors  noted  a  substantial  sensitivity  in  the  degree  of  flow  sep¬ 
aration  and  in  the  unsteady  flow  behavior  to  small  changes  in  the 
constants  of  the  turbulence  model. 
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Chapman  (Ref  39)  summarized,  in  1979,  the  two  key  features  found 
in  successful  solutions  of  the  Reynolds-averaged  Navier-Stokes  solutions 


for  unsteady  turbulent  flows:  (1)  the  flowfield  exhibits  a  single 
dominant  oscillatory  frequency,  (2)  this  frequency  is  an  order  of 
magnitude  lower  than  the  mean  frequency  of  the  turbulent  eddies. 

To  date,  Liou  and  Coakley  have  published  the  only  computational  work 
concerned  with  unsteady  diffuser  flows.  This  differs  from  the  present 
effort  in  that  it  addresses  the  relatively  mild  oscillations  found  in 
the  supercritical  flow  regime.  The  present  effort  is  directed  toward 
the  much  larger  amplitude  instability  in  both  the  internal  and  external 
flowfields  associated  with  the  subcritical  flow  regime. 

1.4  Research  Objectives 

The  review  of  previous  work  is  indicative  of  the  complexity  of  the 
problem.  The  inlet  flow  is  characterized  by  large  regions  of  separated 
flow  and  viscous-inviscid  (shock  boundary  layer)  interaction  in  an 
unsteady,  turbulent  flow.  Experimental  studies  have  been  extremely 
limited  by  instrumentation.  Previous  analytic  works  have  not  fully 
modeled  the  real  flow  and,  as  a  result,  have  arrived  at  incomplete  and 
at  times  contradictory  results.  Today  the  actual  physical  mechanisms 
Involved  in  inlet  buzz  are  still  unresolved. 

The  objective  of  this  effort  is  to  compute  Navier-Stokes  solutions 
for  the  flowfield  about  an  inlet  operating  in  the  near-critical  and 
the  unstable  subcritical  regimes.  These  computations  provide  the  first 
test  of  current  Navier-Stokes  algorithms  to  numerically  capture  the  buzz 
phenomenon  which  occurs  during  subcritical  flow  conditions.  As  such, 
they  offer  the  potential  of  providing  a  new  tool  in  our  effort  to  better 


understand  the  instability  process.  Because  different  instability  mech¬ 
anisms  are  present  in  different  configurations,  one  specific  experiment 
has  been  chosen  for  comparison.  The  external  compression  axi-symmetric 
inlet  tested  by  Nagashima,  Obokata,  and  Asunama  was  selected  as 
the  basis  for  numerical  solutions.  One  steady  state  solution  will  be  ob 
tained  for  verification  of  the  solution  procedure  and  then  several 
solutions  in  the  unsteady,  subcritical  regime  will  be  studied.  Where 
possible,  the  numerical  solution  will  be  compared  with  experiment  and 
analysis  to  provide  a  unified  explaination  of  inlet  buzz. 
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Section  II 


MATHEMATICAL  DESCRIPTION  OF  THE  PROBLEM 


2.1  Governing  Equations 

The  equations  governing  the  motion  of  an  unsteady,  compressible, 
viscous  gas  with  no  body  forces  or  heat  sources  are  derived  from  con¬ 
servation  laws  of  mass,  momentum,  and  energy.  Conservation  of  mass 
(continuity)  states  that  the  rate  of  increase  of  mass  in  a  given 
volume  is  equal  to  the  net  rate  of  inflow  of  mass  through  the  bound¬ 


ing  surface,  or. 


ft  <p> 


(p)  u j  n^  ds 


(2.1) 


making  use  of  cartesian  tensor  notation  with  the  summation  convention. 
Conservation  of  momentum  states  that  the  rate  of  increase  of  momentum 
within  a  control  volume  is  equal  to  the  net  rate  of  inflow  of  momentum 
through  the  surface  plus  the  net  force  acting  on  the  volume  at  the 
surface,  or. 


J ~  (pu±)dv  -  -  J  (pu±)  Uj  n^  ds  +  J x^  n^  ds  (2.2) 

Conservation  of  energy  requires  that  the  rate  of  change  of  energy 
within  a  volume  is  equal  to  the  net  rate  of  transport  of  energy 
through  the  surface  plus  the  rate  of  heat  addition  and  rate  of  work 
done  on  the  volume  at  the  surface. 

J 1^-  (pE)dv  -  -  J (pE)  u j  n^  ds  -  J q^  n^  ds  +  J u±  x^  n^  ds 


with  E  -  ^ei  +  J  \ 

Making  use  of  the  Gauss  divergence  theorem, 


an.  ds 

9  J 


'  1*1 


(a)  dv 


(2.3) 


(2.4) 
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the  equations  may  be  written  in  differential  form,  valid  for  an 
arbitrary  infinitesimal  volume: 

It  (P)  +  ^  <P<9  ■  °.  (2.5) 

h  (oui>  +  jlj-  <pui  uj  -  V  ■  °*  <2-6) 

f?  (oE)  +  <tiuj  E  -  ui  Tij  +  ’j>  ■ 0  (2-7> 

These  equations  are  commonly  known  as  the  Navler-Stokes  equations  and 
describe  the  dynamics  of  very  general  flows.  Complete  derivations 
may  be  found  in  many  books,  for  example,  Liepmann  and  Roshko  (Ref  40: 
332-338) . 

The  stress  tensor  is  composed  of  a  pressure  term  acting  normal  to 
the  surface  and  non-isotropic  term,  the  deviatoric  stress  tensor,  d^, 
producing  normal  and  tangential  stress  as  a  consequence  of  the  viscosity 
of  the  fluid 


ij 


-p6u  +  du 


(2.8) 


For  a  Newtonian  fluid,  the  deviatoric  stress  is  linearly  related  to  the 


rate  of  strain  and  can  be  expressed  as 

8ui.3uj.  ...  3uk 


dij  "  +  +  xsij 

The  heat  flux  vector  is  taken  from  Fourier's  law 


-k 


3T 

ax 


j 


From  Sutherland's  law,  the  viscosity  is 
*--&-)  3/2  (: 


T  +  110°K) 


T  +  110°K) 


X  ■  -2/3  y  (Stoke 's  hypothesis)  and  k  ■  y  C^/Pr 


(2.9) 


(2.10) 


(2.11) 


The  equations  are  closed  by  an  equation  of  state.  For  the  present 
problem,  a  perfect  gas  may  be  assumed: 

p  -  pRT  (2.12) 
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and  if  the  gas  is  calorlcally  perfect. 


ex  -  CvT  (2.13 
Together  with  proper  initial  and  boundary  conditions,  the  equation  set 
(2.5)  to  (2.13)  complete  the  mathematical  problem  description. 

2.2  Axi-symmetric  Coordinate  System 


Since  the  problem  to  be  considered  involves  an  axi-symmetric  inlet 
the  most  natural  coordinate  system  is  a  cylindrical  one.  If  the  inlet 
is  only  examined  at  zero  angle  of  attack  the  flowfield  exhibits  no 
azimuthal  dependence  and  one  of  the  three  momentum  equations  may  be  de¬ 
leted,  with  additional  simplifications  in  the  remaining  equations. 

Equations  (2.5)  to  (2.7)  may  be  generalized  to  any  orthogonal 


curvilinear  coordinate  system  as: 
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For  cylindrical  coordinates 
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and  the  scale  factors  are 
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Substituting,  and  imposing  the  axi-symmetric  requirement 
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the  equations  may  be  expressed  in  vector  form  as: 
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The  vector  components  correspond  to  the  continuity,  two  momentum,  and 
energy  equations  respectively. 


2.3  Reynolds-Averaged  Navler-Stokes  Equations 

In  many  cases,  the  flow  is  turbulent  and  Is  characterized  by  random 
fluctuations  in  all  variables.  Usually,  it  is  impossible  to  numerically 
resolve  all  of  the  length  and  time  scales  involved,  and  so  an  averaging 
process  is  carried  out  to  solve  for  the  mean  motion  which  is  of  primary 
importance  anyway.  According  to  Chapman  (Ref  39)  this  averaging  is 
valid  for  unsteady  flows  as  well,  provided  that  its  time  scale  is  long 
in  comparison  with  the  mean  frequency  of  the  turbulent  eddies,  yet  short 
when  compared  to  the  time  scale  of  the  unsteady  flow. 

Rubesin  and  Rose  (Ref  41)  Introduced  an  averaging  procedure  which 


leads  to  a  convenient  form  of  the  equations.  The  dependent  variables  are 
written  in  terms  of  time  averaged  mean  and  fluctuating  terms: 


and  mass  averaged  mean  and  fluctuating  quantities: 
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The  following  rules  of  averaging  are  assumed  to  apply: 
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In  mass  averaging,  the  mass  flux  is  written  as  a  single  term  and  the 
mean  velocity  is  defined  as 


PU4 
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then 
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which  implies 
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Expressions  similar  to  equations  (2.25)  to  (2.27)  may  also  be  written 
for  the  internal  energy,  e  .  Upon  substituting  the  variables,  in  terms 
of  mean  and  fluctuating  quanlties,  into  equations  (2.5)  to  (2.7)  and 
time  averaging,  the  following  mean  flow  equations  result: 
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The  term 

“i  <Tij  -  2  l>“luP 

has  been  neglected  from  (2.30)  by  order  of  magnitude  considerations. 

The  equations  for  the  stress  tensor  and  heat  flux  vector  assume  the  same 
form  as  their  laminar  counterparts  if  fluctuating  quantities  are  again 
assumed  small  in  comparison  to  mean  values. 

Equations  (2.28)  to  (2.30)  differ  from  their  laminar  counterparts 
only  with  the  exception  of  the  following  additional  terms: 

-  Pu^Uj  (Reynolds  stress)  (2.31) 

puj'ej  (Reynolds  heat  flux)  (2.32) 

These  terms  represent  the  effect  of  the  turbulent  fluctuations  on  the 
mean  motion.  For  example,  the  Reynolds  stress  term  represents  an 
apparent  additional  stress  which  is  a  property  of  the  flow  and  not 


the  medium. 


v,y. 


2.4  Turbulence  Model 


The  turbulent  fluctuating  averages  are  additional  unknowns  for 


which  there  are  no  corresponding  equations.  The  equations  must  be 
closed  at  some  point  by  appropriate  expressions  for  the  unknowns.  Most 
practical  methods  to  date  have  been  based  upon  the  eddy  viscosity  con¬ 


cept.  This  includes  the  simpler  algebraic  models  as  well  as  the  more 
complicated  one  and  two  equation  models.  The  Reynolds  stress  tensor  is 
modeled  as  being  proportional  to  the  laminar  stress  tensor  of  the  mean 


flow,  with  the  coefficient  of  proportionality  being  defined  as  the  eddy 


viscosity 
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A  similar  eddy  dlffusivity  for  heat  flux  can  be  defined  as 
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The  constant  Prandtl  number  is  based  upon  a  Reynolds  analogy,  in  which 
the  eddy  diffusivlties  of  momentum  and  heat  are  related. 

Upon  substituting  expressions  (2.33)  and  (2.34)  into  equations 
(2.28)  through  (2.30)  and  transforming  to  axi-symmetric  coordinates  by 
equations  (2.14)  to  (2.19)  the  Reynolds-averaged  axi-symmetric  Navier- 


Stokes  equations  may  be  written  as: 
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A  Cebeci-Smith  two-layer,  equilibrium  model  (Refs  42,  43,  44)  was 
chosen  to  model  the  eddy  viscosity  in  regions  where  the  flow  was  ex¬ 
pected  to  be  turbulent.  That  is,  in  the  boundary  layer  region  down¬ 
stream  of  the  transition  from  laminar  to  turbulent  flow.  The  eddy 
viscosity  is  assumed  to  be  composed  of  an  inner  and  outer  layer.  In 


both  cases  the  eddy  viscosity  may  be  expressed  as, 
e  -  pil/r(s) 


(2.36) 


where  £  and  1/  are  appropriate  length  and  velocity  scales.  In  the  inner 
region,  Cebeci  (Ref  42)  gives  the  length  scale  for  an  axi-synmetric 
boundary  layer,  with  Van  Driest’ s  modification  as: 
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and  the  velocity  scale  is 


(2.37) 


“  .40  (Von  Kansan  s  constant) 
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y  «  distance  normal  to  the  surface 
yQ  *  body  radius 

u  “  velocity  component  tangent  to  surface 
r(s)  ■  transition  factor 

s  -  distance  along  surface  from  leading  edge 
In  the  outer  region,  the  length  scale  is  proportional  to  the  two  dlmen 
sional  incompressible  displacement  thickness. 

1  "  k2  6i  "  k2 

m  .0168  (Clauser's  constant) 

The  velocity  scale  is  the  reference  velocity 

1/  •  “re£  -  velocity  tar  gent  to  boundary  at  the  edge  of  the 
boundary  layer. 

A  continuous  distribution  of  eddy  viscosity  is  obtained  by  switching 
to  the  outer  expression  when  it  is  first  exceeded  by  the  inner  ex¬ 
pression.  The  equation  used  to  model  the  transition  from  laminar  to 
turbulent  flow  is  taken  from  Dhawan  and  Narashima  (Ref  45) 

T(s)  •  1-exp  (-.412  s2) 

s  -  (s  -  si)/A  si  <_  s  <sf  (2.40) 

A  "  slr-  3/4  ■  s I r  -  i/4 

The  locations  for  turbulence  onset  (s^)  and  fully  turbulent  flow  (s^) 
must  be  known  or  approximated  a  priori. 

The  Cebeci-Smith  model  has  the  Important  advantage  of  its  rela¬ 


tive  economy  compared  with  the  more  complex  differential  equation 


models.  It  has  been  extensively  used  with  good  success  for  attached 
flows  without  large  adverse  pressure  gradients.  Unfortunately  the 
problem  to  be  considered  involves  large  regions  of  separated  flow 
caused  by  highly  adverse  pressure  gradients.  The  author  is  not  aware 
of  any  model  which  is  generally  satisfactory  for  unsteady  separated 
flows.  The  algebraic  model  assumes  an  equilibrium  between  production 
and  dissipation  of  the  turbulence  which  is  not  true  in  regions  of 
rapid  change  such  as  adverse  pressure  gradients  and  separated  boundary 
layers.  Relaxation  and  pressure  gradient  corrections  have  been  used 
in  these  regions  to  improve  agreement  with  experiment  but  usually 
require  model  adjustment  for  each  flow  field.  Another  common  cor¬ 
rection,  intermittency ,  was  not  implemented  due  to  its  little  influence 
on  the  flow  field  and  the  eventual  merger  of  the  downstream  boundary 
layers.  Since  transition  must  be  modeled,  the  location  and  extent 
of  the  transition  region  should  be  known  from  experiment  or  otherwise 
rationally  estimated.  Finally,  what  part  of  the  turbulent  spectrum 
should  be  modeled  when  the  subject  of  the  calculation  is  itself  a 
periodic  instability  must  be  considered.  In  essence,  this  means  that 
the  limits  over  which  the  turbulent  fluctuations  are  averaged.  Equation 
(2.23),  should  be  appropriate  for  the  problem  at  hand.  If  interest 
is  focused  on  the  direct  calculation  of  the  flow  instability  and  this 
is  resolvable  to  a  limited  extent  on  the  computational  grid,  the 
period  over  which  averaging  is  performed  and  the  eddy  viscosity  co¬ 
efficient  which  models  the  fluctuating  averages  should  correspond  only 
to  the  higher  frequencies  which  are  not  resolvable  (Ref  46) .  The 
Cebeci-Smith  model  was  adopted  as  a  baseline  to  provide  an  approximate 
representation  of  the  turbulence  at  a  reasonable  cost. 


Section  III 


NUMERICAL  METHODS 


3.1  Body  Oriented  Coordinate  System 

Coordinate  Transformation  to  Computational  Plane 
The  governing  partial  differential  equations  must  be  discretized 
over  a  finite  number  of  points  at  which  approximate  solutions  are  gen¬ 
erated.  A  body-oriented  coordinate  system  is  usually  chosen  to  align 
the  coordinate  axes  with  the  body  contours  in  order  to  simplify  the 
implementation  of  boundary  conditions.  Such  a  grid  also  provides  suit¬ 
able  resolution  in  regions  of  rapid  change  through  coordinate  grid 
stretching.  The  solution  algorithm,  however,  must  be  on  a  uniform  grid 
if  standard  finite  difference  formulas  are  to  be  valid.  This  is  accom¬ 
plished  by  transforming  the  physical  coordinates  (x,y)  to  computational 
coordinates  (£ , n) • 

If  the  functions  5(x»y)  and  n(x,y)  and  their  first  partials  are 
continuous  and  the  Jacobian,  J,  4  0,  then  the  coordinate  transformation 
is  admissible.  The  following  relations  may  then  be  obtained  by  applica¬ 
tion  of  the  chain  rule: 
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The  two  Jacobian  matrices  are  inverses  and  thus  the  transformation 


derivatives  may  be  evaluated  as: 
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Equation  (2.35)  may  be  expanded,  using  equation  (3.1)  to  give 
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The  equations  are  usually  written  in  divergence,  or  conservation 
law,  form,  (Refs  47,  48,  49,  50,  51)  because  of  its  ability  to  capture 
"weak  solutions."  This  is  a  consequence  of  the  exact  cancellation  of 
intercell  fluxes,  which  assures  strict  conservation  of  mass,  momentum, 
and  energy.  Conservation  form  is  most  advantageous  for  flows  involving 
shocks.  This  form  can  be  shown  to  satisfy  the  Ranklne-Hugoniot  shock 
jump  relations,  ensuring  correct  shock  location,  strength,  and  speed. 


Upon  multiplying  by 


J"1  y 


Equation  (3.5)  may  be  rewritten  as: 
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It  is  not  difficult  to  show  that  the  terms  on  the  right  hand  side  are 

both  identically  zero.  The  remaining  source  term  is  a  consequence  of 

the  ignorable  coordinate,  @,  and  must  remain  in  the  axi-symmetric 

equations.  If  the  contravariant  velocity  components 
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are  defined,  the  governing  equations  may  be  written  in  a  compact  form 
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The  superscripts  for  averaged  variables  have  been  dropped  and  par- 
tials  with  respect  to  x  and  y  in  the  shear  stress  terms  are  evaluated  by 
application  of  Equation  (3.1).  For  example, 

Txy  *  <vi+£X-^  5y  +  ^  ny  +  35  +  lr\  V  (3*9) 

This  is  the  form  of  the  equations  used  for  numerical  integration. 


Coordinate  Generation 

Two  different  approaches  were  used  to  generate  coordinate  systems 
for  the  axi-symmetric  inlet  considered.  A  popular  approach  is  to 
map  a  coordinate  system  about  any  arbitrary  body  such  that  resulting 
computational  plane  is  then  a  unit  square.  This  method  has  the  out¬ 
standing  advantage  that  the  solution  algorithm  is  independent  of  any 
particular  geometry.  Since  both  the  internal  and  external  flow  about 
the  inlet  must  be  computed,  the  grid  must  be  wrapped  about  the  cowl 
lip.  The  resulting  C-shaped  grid  is  shown  in  Figure  2.  Thompson's 
method  (Ref  52)  was  used  to  generate  the  grid  as  a  solution  to  the 


elliptic  system  of  equations; 

v25-4+4-° 
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The  forcing  function,  T(5,n)»  which  determines  coordinate  line  attraction 
near  the  boundaries  was  given  by  Knight  (Ref  53)  and  iteratively  correct¬ 
ed  to  maintain  a  constant  step  size  at  the  boundaries.  This  approach  was 
used  to  generate  acceptable  flowfield  solutions  for  the  upstream  inlet. 

It  was  abandoned,  however,  because  of  several  distinct  disadvantages. 

The  cowl  lip  tested  experimentally  was  sharp.  To  avoid  the  singularity 
in  the  transformation,  the  lip  must  be  rounded.  As  the  radius  of  cur¬ 
vature  becomes  smaller,  the  transformation  derivatives  become  larger  anu 
the  allowable  time  step  in  the  solution  algorithm  becomes  smaller.  Any 
finite  radius  of  curvature  also  produces  a  detached  lip  shock  which  is 
not  present  with  a  sharp  lip  at  the  condition  tested  experimentally. 

Also  since  the  total  number  of  n  grid  lines  is  the  same  in  both  the 
interior  and  exterior  regions  of  the  inlet,  a  trade  off  must  be  made 
between  the  coarse  grid  exterior  to  the  inlet  with  a  larger  truncation 
error  and  a  very  dense  interior  grid  resulting  in  a  lower  allowable 
time  step  in  the  solution  algorithm.  Finally,  it  is  difficult  to  de¬ 
termine  forcing  functions  which  counter  the  natural  tendency  of  the 
elliptic  system  to  cluster  points  about  convex  surfaces  and  away  from 
concave  surfaces. 

A  much  better  approach  is  to  map  the  physical  plane  into  an  L-shaped 
domain,  Figure  3,  in  which  the  cowl  lip  is  treated  as  an  interior  region. 
The  cowl  tip  is  resolved  to  within  one  grid  cell  and  the  lip  need  not  be 
rounded.  Such  a  mesh  may  be  quickly  and  inexpensively  produced.  An 


Physical  Plane 


If 


Computational  Plane 


CV*. 


Fig.  2  C-Grid  Coordinate  Systea  -  Inlet  Forebody 
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Fig.  3  L-Grid  Coordinate  System  -  Inlet  Forebody 


algebraic  mesh  generation  scheme  was  used  which  clustered  points  near  the 
centerbody  and  cowl  walls  by  an  exponential  stretching,  smoothly  blended 
to  a  constant  step  size  in  the  center  of  the  duct.  The  external  region 
of  the  grid  was  generated  between  the  inner  grid  boundary  and  the  specified 
outer  boundary  with  a  region  of  exponential  stretching  blended  to  a  con¬ 
stant  step  size  region.  First  derivative  continuity  was  maintained  through 
out  the  grid.  The  grids  generated  for  the  full  inlet  configurations  are 
shown  in  Section  IV. 

Once  the  grid  has  been  generated,  the  transformation  derivatives 
are  computed  and  stored  once  and  for  all  using  Equations  (3.3).  Stand¬ 
ard  second  order  accurate  difference  formulas  are  used  to  obtain  de¬ 
rivative  values.  In  the  interior,  central  differences  are  used 


l+l.J  ‘1-1, 

2A£ 


,1  ill  .,2. 

35  i.i 


"  fl,J-l  .  (i  i!i 
2A*  an3. 


(3.11) 


where  f  ■  x,y  and  the  leading  truncation  error  term  is  indicated  in  paren¬ 
thesis.  On  the  boundaries,  including  the  cowl  lip,  transformation  deriva¬ 


tives  are. similarly  given  by  one  sided  forward. 


and  backward  differences 


3.2  Solution  Algorithm 

MacCormack's  Method 

MacCormack's  explicit  finite  difference  algorithm  (Refs  54,55)  was 
chosen  for  numerical  Integration  of  the  Navier-Stokes  equations  (3.8), 
based  upon  its  proven  success  in  a  wide  variety  of  problems.  The  method 
provides  consistent  time  accurate  resolution  for  unsteady  flows  and  is 
noted  for  its  robustness  and  excellent  shock  capturing  ability. 

The  method  is  best  illustrated  by  consideration  of  a  simple  exam¬ 
ple  in  non-transformed  variables.  Since  multi-dimensional  systems  of 
equations  may  be  split  into  sequences  of  one-dimensional  operators  by 
the  method  of  fractional  steps,  the  scheme  may  be  demonstrated  by  con¬ 
sidering  the  non-linear  set  of  equations 


^(u)+|j(r).o 


(3.14) 


where  F  ■  F(U)  and  U  *  U(t,x)  are,  in  general,  vectors. 


The  equation  may  be  discretized  in  time  and  space  as 


F"  -  F  (U  (nAt ,  iAx) ) 


x  ■  iAx  ,  t  “  nAt  , 

A  second  order  time  accurate  predictor-corrector  algorithm  may  be 
written  as 


U 


rn+l 


„n  . .  .  _n 
U.  -  At  3  F. 
i  x  i 


(3.15) 


U"+1  -  u"  -  \  AtOxF"  +  3X^+1>  (3.16) 

which  is  the  modified  Euler  predictor-corrector  method.  The  notation, 
n+1,  is  used  to  denote  a  provisional  n+1  time  level  given  by  the  predict 
or.  MacCormack's  algorithm  results  when  the  spatial  derivatives  are 


alternately  given  by  one  sided  first  order  differences 


uj*1  -  Uj  -  At  7x  Fj 


uf  1  -  -  |  *t  (Vx  F"  +  Ax  Ff1) 


(3.17) 

(3.18) 
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Upon  substituting  from  Equation  (3.17),  Equation  (3.18)  may  be  written  as 


of1  -  ±  [«?  +  uf 1  -  it  ix  F"Ti] 


(3.19) 


The  result  in  a  second  order  accurate  in  space  and  time,  conditionally 
stable,  procedure  of  Lax-Wendroff  type  with  effectively  central  spatial 
differences.  This  is  shown  in  Appendix  A,  which  summarizes  the  diff¬ 
erence  operators  used  in  MacCormack's  algorithm.  The  scheme  is  dissi¬ 
pative  due  to  the  asymmetry  in  the  differencing  and  this  contributes 
much  to  the  robustness  of  the  method.  Appendix  B  demonstrates  rig¬ 
orously  the  second  order  accuracy  of  the  scheme  in  which  the  difference 
equation  generates  an  approximation  to  the  differential  eqn  of  the  form 


ut  +  Fx  +  0( At2)  +  0  (AtAx)  +  0  (Ax2)  -  0 


(3.20) 


The  Navier-Stokes  equations  contain  second  derivatives  resulting 
from  the  viscous  stress  terms.  These  appear  in  the  F  and  G  terms  as 
first  derivatives.  If  the  outer  differencing  of  the  F  and  G  terms  is 
forward  (backward)  and  the  interior  derivatives  occur  in  the  same  vari¬ 
able  as  the  outer  derivative,  these  terms  are  then  evaluated  by  first 
order  backward  (forward)  differences.  Mixed  derivative  terms  are 
evaluated  by  2nd  order  central  differences  for  the  interior  terms. 

As  an  example,  a  component  of  F  having  the  form 
,3u  ,  3v. 

F  " 

is  differenced  by  Equation  (3.16)  as 

u"+J  -  U®  +  At  V  Fj  -  u"  -  ^  (Fj  -  F“  .  )  (3.21) 

i,j  i.j  x  i,  j  i,j  Ax  i,J  i-l.J 


„n  n  /  (i  n  .  .  n  v 

Fl ,3  '  “i.J  <Sy  “1,3  +  4*  vi,3) 


_n  n  s  ^  q  .  .  a  \ 

Fi-l,l  ■  “i-l.j  <5y  “1-1,3  +  Vl,3) 


(3.22) 


and  by  Equation  (3.19)  as 


t. 


.V-.  ’ 


pps 


;«r» 

*  Cj 


<5  ■  I  +  -  “  4x  O 

- 1  "ij +  -  f  -  O 


(3.23) 


with 


r.n+1 


n+1 


n+1 


— ■  A  u  •  A  /  n  IX  I  X  1.7  n+1  v 

Fi+l,j  "  Vi+l,j  (6y  ui+l,j  7x  Vi+l,j) 


(3.24) 


r.n+1 


nH  *  *  n+1  *  p  n+1  ,  _  n+1* 

Fi.J  "i.j  (5y  ui,j  +  7»  vi.3) 


The  result  Is,  again,  an  effectively  central,  second  order  accurate 


difference  approximation  to  the  second  derivative  terms. 


2  2 
3  u  .  3  v 

”d^T 


If  viscosity  is  considered  constant.  Equation  (3.23)  may  be  rewritten, 


using  the  notation  of  Appendix  A,  as 


,,n+l  _  „n  At  .  /t7  n  .  .  .n+1* 

Vj  "i.j  •  "T  V’x  ui,j  +  4x  "i.j1 


At  _  .  ,  n  n+1 \ 

-  v~2  Vx  <Vl,j  +  *i,J> 


(3.25) 


in  which  the  structure  of  the  differencing  is  transparent. 


Splitting 

For  multi-dimensional  systems  of  equations,  additional  efficiency 
may  be  gained  if  the  equations  are  split  into  a  sequence  of  one  dimen¬ 
sional  operators.  The  two  dimensional  equation 


h  <“>  +  It  (f)  +  k (G)  ■ 0 


(3.26) 

may  be  advanced  in  time,  using  the  method  of  fractional  steps  (Refs  54, 
55,  56,  57,  58)  by  the  symmetric  operator  sequence: 

(3.27) 


1,0+1  "  Lv  <%”/2  L  (At)  L(%“/2  Un 
y  ©  x  y  © 
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1  „  *  **  ** 

2  «i.J  +  B1.J  -  “  \  F1,J> 


Splitting  is  efficient  when,  because  of  stability,  the  allowable  time 

step  in  one  direction  is  much  smaller  than  the  other.  That  is,  where  m 

the  largest  integer  ratio  of  allowable  time  steps  in  the  L  and  L 

x  y 

operators,  is  greater  than  two.  The  solution  may  be  advanced  several 
times  in  one  operator  before  the  other  need  be  advanced.  This  sit¬ 
uation  typically  occurs  in  the  boundary  layer  region,  where  in  order 
to  resolve  the  large  flow  gradients  normal  to  the  wall 
Ay/ Ax  <<  1 

Lomax  (Ref  58)  has  shown  that  if  a  symnetric  operator  sequence  is 
followed,  the  split  algorithm  retains  full  second  order  temporal  and 
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spatial  accuracy.  This  Is  outlined  In  Appendix  C 


Computational  Coordinates 


MacCormacks  method,  as  outlined,  may  now  be  applied  directly  to  the 
transformed  Navier-Stokes  Equations  (3.8).  The  solution  is  then 


advanced  in  time  by  the  operator  sequence. 

=  L  (.At  )  L.(At_)  L  (At)  /  . 
i,j  n  n  £  £  n  n  i,j 


(3.32 


where  m  -  -—*■  =  2,  £  ■  (i-l)A£,  n  *  (j-l)An,  t  «  nAt 
n  5 

and  i  **  1,2,3,.  ..,ik  j  -  1,2,3,  ...,jk 


The  L  operator , 

n 

U*  «  L  (At  )  Uj 
i»j  n  h  i,j 

updates  the  radial  split  equation 


(3.33 


It  (U)  +  h  (G>  ■  H 


(3.34 


X  Ail- 

U .  4  4  -  T11  4  -  G"  ,  -  )  +  At  H” 

i, j  i.J  An  i,j  i,j-r  n  i,j 


7  .  ■  4  [U°  .  +  U*  •-  7— 11  (G*  ...  -  G"  )  +  At  h"  ]  (3.35 

i,j  2  L  i,j  i , J  An  i, j+1  i,j  n  i,j 


*  * 


The  Lj.  operator, 

u‘i!j  ■  w  uIj 

updates  the  axial  split  equation 


(3.36 


k(U)+k (F) ' 0 


(3.37 


U**  m  / - -S-  (F*  -  F*  ) 

Ui,j  Ui,j  AC  vri,j  ri-l,J; 


U**.  -  i  [U*  +  U**.  -  ^5.  <F**.  . 


-  F**.)l 


(3.38 


The  method  carries  through  In  a  straight-forward  manner,  with  the  dif¬ 


ferencing  now  being  either  forward,  backward,  or  central  in  the  computa¬ 


tional  variables  £,q.  However,  the  interior  derivatives  in  the  source 


term,  tf,  are  approximated  with  central  differences  in  both  £  and  n  for 


both  the  predictor  and  corrector  in  the  radial  equation. 


Stability 


MacCormack’s  method  is  conditionally  stable  with  a  corresponding  re¬ 


striction  on  the  allowable  time  step.  Stability  of  a  given  algorithm  is 


difficult  to  ascertain.  The  most  successful  approach  to  date  is  that  of 


Von  Neumann  (Refs  50,51).  The  governing  equations  must  be  linearized  and 


the  coefficient  matrices  assumed  to  be  locally  constant  or  slowly  vary¬ 


ing.  Consistent  with  the  linearized  form  of  the  equations,  the  growth  or 


decay  of  each  Fourier  component  is  examined  independently.  The  approach 


is  to  assume  a  solution  of  the  form 


v"  ,  -  ,”<kr,  k  )e  I(k5145  +  V4n) 

*■»  3  s  n 


(3.39) 


where 


I  -  /  -1 


<p  (k^,  k^)  -  amplitude  function  at  nth  time  level 


k^  -  wave  number  in  £ 


k  -  wave  number  in  n 

n 


A  difference  equation  may  be  obtained  by  applying  the  two  step  MacCormack 


predictor-corrector  to  the  linearized  equations.  Provisional  (  ) 


values  may  be  eliminated  in  the  corrector  by  substitution  from  the 


predictor.  Equation  (3.39)  is  then  substituted  into  the  resulting  cor¬ 


rector  equation.  Considerable  manipulation  will  lead  to  an  equation  of 


the  form 


kv)  p  G(kg ,  y  *n(y  y 


(3.40) 


from  which 


$n(k?,  \)  “  [G(k5,  kn)  ]n  <^°(kc,  kn)  .  (3.41) 

For  a  fixed  time,  t  =  nAt,  as  the  time  step  is  refined,  the  iteration 

number  approaches  the  limit 

Lim  (|^)  *  » 

At  -*■  0 

Thus,  for  stability,  the  matrix  norm  of  the  amplification  matrix,  G,  must 
remain  bounded, 

I  |  Gn  |  |  <  M  as  n  -*■  ®. 

It  can  be  shown  that  a  necessary  condition  is  that 

|X|  <_  1+0  (At),  (3.42) 

where  A  is  the  largest  eigenvalue  of  the  amplification  matrix  for  all 


possible  combinations  of  wave  number  k^ ,  k^.  This  is  frequently  termed 
the  spectral  radius. 

Because  the  equations  are  solved  in  factored  or  split  form,  the 
stability  of  the  L ^  and  L ^  operators  may  be  ascertained  Independently. 
Upon  rewriting  the  non-conservative  form  of  the  governing  equations  (3.5) 
in  the  more  convenient  variables 


(3.43) 


and  linearizing,  the  axial  split  equation  becomes 

£  +  +  V° 


(3.44) 


and  the  radial  is 

£♦»£♦»£  0.45, 

an 

A,  B,  C,  D,  E^,  Ej,  F^,  F ^  are  4x4  constant  matrices  which  are  coefficients 


of  the  inviscid,  diffusive,  mixed  derivative  and  source  terms  respec¬ 
tively.  FjV  and  F^V  appear  because  of  the  axi-symmetric  form  of  the 
equations  and  do  not  appear  in  the  two  dimensional  case.  Richtmyer  and 
Morton  (Ref  50:  Sec  8.4)  demonstrate  that  stability  is  unaffected  by 
the  source  or  lower  order  terms  so  they  may  be  safely  Ignored.  Con¬ 
sequently,  stability  conditions  for  the  axisymmetric  equations  are  the 
same  as  thd^two  dimensional  stability  bounds. 

MacCormack  and  Baldwin  (Ref  55)  have  analyzed  separately  the  stability 
of  the  inviscid,  diffusive  and  mixed  derivative  parts  of  Equations 
(3.44)  and  (3.45)  by  the  von  Neumann  method  for  the  non-transformed 
equations.  The  dominant  stability  component  is  the  Inviscid  limit,  with 
smaller  viscous  corrections. 


The  algebra  involved  is  quite  lengthy  and  need  not  be  detailed  here, 
is  of  interest  are  corresponding  stability  restrictions  for  the  trans 
formed  equations.  Since  the  transformed  equations  are  identical  in 
form  with  those  considered  by  MacCormack,  Appendix  D  demonstrates  how 
similar  stability  restrictions  may  be  generalized  for  the  transformed 
equations.  The  maximum  allowable  time  step  is  then 

At_  «  2  a  min  [min  (At  ),  4  min  (At-  )] 

6  M  ni,J  i.J 


where 


At  min  [min  (At  ),  4  min  (At-  )] 

n  i»j  ni,j  i.J  j 


a  ■  0.  -*■  1.0  (usually  .9) 


|Ur|  +  c  +[ 


AS 


At 


ni.j  20  €L 

l\i  +<:  +  tir  +  ir-]/P 

n  n 


i.j 


Also, 


AS 


5  r 


AL¬ 


LS 


An 


'5  2+C  2 
x  y 


/  2  2 

\  +r>y 


U, 


U 

c 


'e  2+s  2 

x  y 


0X  -  max  [  1 2  (y+e)  +  (X+Xt)  | ,  y(|^  +  §^-)] 
02  »  »;|(y+e)(X+XT)| 


(3.49) 


(3.50) 


This  same  stability  condition  was  given  by  Knight  (Ref  53)  in  his  solu¬ 
tion  for  the  two-dimensional  transformed  equations.  For  many  viscous  flows 
the  time-step  restriction  in  the  operator  is  quite  severe.  The  equa¬ 
tions  become  very  stiff  in  the  sense  that  the  allowable  time  step  is  much 
smaller  than  what  is  required  to  accurately  resolve  the  temporal  dynamics 
of  the  flow.  It  is  the  desire  to  avoid  this  restriction  that  motivates 
the  use  of  an  implicit  method  in  which  spatial  derivatives  are  evaluated  at 
the  n+1  time  level  with  unconditional  stability.  References  (59),  and  (60) 
are  representative  of  current  implicit  methods.  The  basic  explicit  method 
was  used  for  the  present  work,  however,  because  of  its  robustness,  its  time 
accurate  resolution,  and  its  easy  implementation  on  a  vector  processor. 


Numerical  Damping 

MacCormack's  algorithm  is  a  ''shock-capturing"  technique  in  which  a 
discontinuity  such  as  a  shock  must  be  resolved  over  several  mesh  points 
rather  than  the  order  of  a  mean  free  path  length  in  which  it  physically 
occurs.  A  solution  in  this  region  will  normally  exhibit  non-physical 
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oscillations  before  and  after  the  shock.  At  best,  these  are  unappealing 
and  they  nay  lead  to  numerical  Instability.  MacCoraack  (Ref  55)  has 
Introduced  a  fourth-order  smoothing  term  which  is  effective  in  damping 
the  shock  induced  oscillations  and  negligible  elsewhere.  Damping  is 


incorporated  into  the  L ^  operator,  equation  (3.38)  as 


i»j 


_*  * 

Fi.j +  vi.i 


_*  *  * 
Fl-l,j  “*  Fl-l,j  + 


**  **  ** 
Fi+1  fi+l,j  +  yi+l,j 


(3.51) 


**  **  ** 

Ft,r  + 
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+*/5xi,j+5yi,j  ci,j 
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x(ui+l,j"tii,j) 


1+1,  j 
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'1+1 
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,j  V  xi+i,j  yi,: 
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1+1,  J 


U**  ) 
i.j' 


(3.52) 


8  *  0.0  -*■  5.0  (usually  1.0  -*•  2.0) 

The  net  effect  is  an  eddy  viscosity  term  of  the  form 


,  *  U  +  r'  2,-2 
„  .r3  3  r  c  C+C  c 

8  At  LK  Jjr  L  _  X  y 


4p 


£jl 
35 2 


3(1  "j 
35  J 


(3.53) 


added  to  the  difference  Equation  (3.38).  Damping  is  similarly  incorpor¬ 
ated  into  the  L  operator.  Equation  (3.35). 


3.3  Boundary  and  Initial  Condition  Implementation 


Boundary  Conditions 

Boundary  conditions  must  be  updated  after  each  pass  through  the 
predictor  and  corrector  for  each  operator.  Referring  to  Figure  (3), 
they  may  be  classified  in  several  different  categories: 

Solid  wall  -  A,  C,  (and  B  for  no  mass  flow  through  the  nozzle) 
Inflow  -  F  , 

Outflow  -  B,  D  , 

Outer  Boundary  -  E  , 

Centerline  -  A,  ahead  of  centerbody. 


Solid  Wall  Boundary  Conditions 


The  velocity  components  are  set  by  the  no-slip  condition  as 
u  *  v  *  0  . 


(3.54) 


Temperature  is  assumed  to  be  the  adiabatic  wall  temperature  specified 


from  the  free  stream  condition  as 


T  -  (T  ).  [l-r(s)]  +  (T  )„_  T(s) 
w  aw  L  J  aw  TB 


(3.55) 


where 


<Ta„>I.  *  T-  +  >¥?  2^ 


(I,»>TB  ■  T-  +  ^  2C, 


A  compatibility  relation  is  used  to  update  the  wall  static  pressure  by 
finding  the  normal  component  of  the  two  momentum  equations  at  the  wall. 
After  applying  the  no-slip  condition  and  neglecting  the  shear  stress  terms 
by  order  of  magnitude  considerations,  the  compatlblity  condition  becomes 


i£  .  r(iE  r  +  iE  n  )  £  +  (i£  r  +  iE  n  }  -T]  .  r 
3n  U3E  Sc  3n  V  S  3n  nv;  L 


\i  +  \  1. 


r~2  2 

\  +  ny 


(3.56) 


O'. 


This  may  be  discretized  (on  the  centerbody,  for  example)  as 


A  P.  . 
n  1,1 


C  n  +5  n 
x  x  y  y 
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n  +  n 
x  y 
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Sc  pi,2 


(3.57) 


to  update  the  wall  pressure  from  the  previously  computed  values  at 
neighboring  interior  points.  Once  u,  v,  T,  and  p  are  known,  all  other 
variables  may  be  computed  from  these  known  values. 


Inflow  Conditions 

The  inflow  boundary  is  specified  to  be  the  unperturbed  freestream 
condition: 

u 
v 
T 
P 

» 

Since  the  flow  is  supersonic,  there  is  no  downstream  influence  on  the 
inflow  boundary.  Some  care  mist  be  taken  In  the  placement  of  the  upstream 
boundary  to  ensure  the  freestream  condition  is  valid.  In  particular,  the 
upstream  boundary  should  be  well  ahead  of  the  most  forward  travel  of  the 
normal  shock  during  the  buss  cycle  or  the  shock  will  reflect  off  the 
boundary . 

Outflow  Boundary  Conditions 

The  no-change  or  seroth  order  extrapolation  was  used  to  specify  the 
outflow  based  upon  the  computed  Interior  values  as  follows 


u 

00 

o 

T 

00 

p. 


(3.58) 


(3.59) 


Using  1st  order  one  sided  differences 


The  flow  at  the  downstream  boundaries  Is  also  supersonic  except  in  the 
boundary  layer  regions.  As  a  consequence  of  the  domain  of  dependence  in 
supersonic  flow,  the  downstream  boundary  has  no  effect  upon  the  interior 
of  the  flow. 

The  Outer  Boundary 

The  outer  radial  boundary  was  placed  such  that  the  shock  would  exit 
through  the  outflow  boundary  in  both  the  steady  state  and  buzz  cases.  The 
outer  boundary  would  then  be  set  as  the  unperturbed  freestream  value. 
However,  during  its  forward  travel,  the  strong  normal  shock  could  briefly 
impinge  on  the  upper  boundary  and  reflect  into  the  domain.  To  avoid  this 
problem,  a  no-reflection  condition  was  used  which  allowed  the  shock  to  exit 
without  reflection  and  otherwise  recovered  the  freestream  value. 

Flow  conditions  on  the  outer  boundary  are  extrapolated  from  the  in¬ 
terior  along  the  outward  running  characteristic  as 

(3.61) 

This  equation  is  strictly  valid  only  for  two  dimensional,  inviscid, 
supersonic,  steady,  isentroplc  flow  in  a  simple  wave  region  in  which 
there  is  no  change  in  the  inward-running  characteristic.  Along  the  out- 
ward-running  characteristic  then,  the  flow  properties  remain  constant. 

This  relation  will  also  be  valid  for  axi-symmetric  flow  at  large  radius 
and  additionally  has  proven  to  be  a  reliable  boundary  condition  for  more 
general  flows  as  shown  by  Roache  (Ref  51:  282) . 


Referring  to  Figure  4,  the  condition  is  implemented  as  follows: 


(1.)  compute  the  outward  running  characteristic  orientation  at 
the  grid  point  (l,jk),  from  the  flow  angle  and  Mach  angle 

8  ■  tan’1(o)i,ik 


a  ■  tan 


2  ,  'i,Jk 


M-l 

(2.)  extrapolate  backward  to  find  the  intersection  with  the  line 
n(Jk-l)  ■  const  as  a  solution  of  the  equations 

yc  “  y(L»  Jk-1)  -  (x(L,  jk-1)  -  xc) 


yc  ■  y(i,jk)  -  tan  (a+0)  (x(i,jk)  -  x£) 


(3.)  linearly  interpolate  along  arc  length  between  known  Interior 
grid  points  along  n(jk-l)  Line, 

V(i,jk)  -  V(xc,  yc)  -  V (L-l ,  jk-1)  (V(L, jk-1)  -  V(L-l,jk-l)) 

(3.62) 


where 


[xc-x(L-i,j-i)r  +  Cyc-ya-i,j-i)]^ 


As  -  /.  2  2 

Ax  +Ay 


Fig.  4  Characteristic  Boundary  Condition  Implementation 
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Centerline  Boundary  Conditions 


s'. 


! 

r. 

r. 

i*. 


For  an  axi-symmetric  flow,  along  the  centerline  a  symmetry  condition 
must  exist: 


(3.63) 


and 


v  -  0 

This  is  implemented  as  a  one  sided  difference  in  the  transformed 
variables  as, 

A  VT  .  -  0  (3.64) 

n  i,i 


Initial  Conditions 

The  solution  procedure  requires  that  some  initial  flow  field  be 
specified.  This  is  then  integrated  foward  in  time  to  a  steady  state 
which,  if  it  exists,  is  Independent  of  whatever  initial  conditions  are 
specified.  The  situation  is  less  clear  for  the  unsteady  cases  considered. 
A  purely  periodic  solution  would  be  expected  to  be  independent  of 
initial  conditions  as  well.  The  buzz  phenomenon  is  quasi-periodic  with 
a  definite  period  produced  by  the  dominant  lower  frequencies  but  not 
repeatable  due  to  the  superposition  of  non-commensurable  higher  frequency 
components.  A  solution  in  this  case  is  not  necessarily  independent  of 
the  specified  initial  conditions.  This  is  especially  true  when  only 
a  few  cycles  are  computed  and  the  flow  has  several  possible  modes. 

Several  different  types  of  initial  conditions  were  used  for  the 
cases  presented  in  Section  IV.  These  are  briefly  outlined  below. 


Near-Critical  Steady 

Upstream  flow  conditions  (from  inflow  boundary  to  some  distance 
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downstream  of  cowl  lip)  were  specified  as  uniform  freestream  flow.  Down¬ 
stream  conditions  were  taken  from  one  dimensional  flow  relationships 
based  upon  area  ratio,  with  an  assumed  location  of  an  ideal  normal  shock 
in  the  diffuser. 

Unstable  Subcrltical  With  Hass  Through-Flow 
Here  the  converged  steady  state  solution  was  used  as  initial  condi¬ 
tion  for  the  subcrltical  flow  with  an  impulsive  decrease  in  the  downstream 
throat  area. 


Unstable  Subcrltical  With  No  Mass  Through-Flow 
Upstream  flow  conditions  were  taken  from  the  previous  unsteady  results 
in  which  the  normal  shock  has  been  forced  from  the  diffuser.  Within  the 
diffuser,  a  no  flow  condition  was  Imposed  (u,v-0)  Static  pressure  was 
fixed  at  the  experimentally  determined  point  of  maximum  pressure  recovery 
before  the  onset  of  buzz.  Temperature  was  also  fixed  uniformly  at  the 
adiabatic  wall  temperature,  which  is,  of  course,  close  to  the  stagnation 
temperature. 

3.4  Implementation  on  a  Vector-Processor  (Cray  1-S) 

All  calculations  presented  in  this  report  were  done  on  the  Cray  1-S 
at  the  Air  Force  Weapons  Laboratory,  Kirtland  Air  Force  Base,  Albuquerque 
New  Mexico.  Access  was  through  the  Scientific  and  Management  Network 
(SAMNET)  subsystem  of  the  ARPA  Network  (ARPANET)  and  also  by  long  dis¬ 
tance  commercial  lines  dial  up. 

The  Cray  1-S  is  known  as  a  vector  processor  and  uses  a  pipeline 
architecture  to  dramatically  increase  the  execution  rate  of  a  given 
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program.  An  interesting  non-technical  introduction  to  the  subject  of 
vector  processors  and  "super  computers"  has  recently  been  given  by  Levine 
(Ref  61) .  In  essence,  a  pipeline  architecture  continously  feeds  new 
operands  into  the  ’pipeline"  in  which  subsequent  elements  of  a  vector 
operation  immediately  follow  the  initial  element  as  it  is  broken  down 
into  segments  and  proceeds  to  completion.  After  the  time  required  for 
the  initial  operand  to  pass  through  the  pipeline,  T8tartup*  a  new 
result  will  be  generated  every  clock  period,  Tq,  thereafter.  A  vector 
operation  of  length  N  will  then  execute  in  the  time,  (Ref  62), 


T  -  T  +NT 

tot  startup  o 


(3.65) 


The  Cray  1-S  has  a  peak  execution  rate  in  excess  of  100  million  floating 
point  operations  per  second  (MFLOPS) . 

Vectorization  is  automatically  done  by  the  FORTRAN  compiler  (Ref 
63:  Chapter  4)  if  the  structure  of  the  code  is  such  that  vectorization 
is  not  inhibited.  The  Cray  compiler  vectorizes  the  innermost  do  loop 
in  which  the  vector  length  is  equal  to  the  difference  between  starting 
and  stopping  values  of  that  loop  index.  Vectorization  is  inhibited  by: 

(a)  Recursive  relationships  in  the  inner  loop  index  , 

(b)  Subscripts  which  are  not  constant  integer  increments 
of  the  loop  index  , 

(c)  60  TO,  IF,  Call  or  other  logical  call  statements  , 
occurlng  in  the  innermost  loop.  An  efficient  vectorized  code  requires 
vectorizable  do  loops  with  a  long  vector  length  so  that  execution 

time  is  not  dominated  by  the  startup  time.  The  execution  rate  approach¬ 
es  an  asymptotic  limit  near  integer  multiples  of  the  vector  register 
length  of  sixty  four. 


MacCormack's  explicit  algorithm  is  easily  vectorized  and  this  has 


been  done  for  the  present  code.  Vector  length  was  set  by  the  1  variable 
corresponding  to  the  £  coordinate.  The  minimum  vector  length  is  then 
near  fifty  and  the  maximum  over  two  hundred.  A  common  measure  of  the 
efficiency  of  a  given  algorithm  on  a  given  computer  is  the  rate  of  data 
processing 


RDP 


(CPU  Time) 


(AGrid  Points) (Alterations)  . 

In  the  present  case,  a  data  processing  rate  per  operator  of 
(RDP) 


(3.66) 


1.166x10  ^  sec 


oper 

has  been  achieved.  This  is  in  good  agreement  with  the  results  reported 
by  other  users  of  vectorized  explicit  codes  on  the  Cray  1-S,  and  Cyber 
203,  (Refs  28,64). 

It  should  be  noted  that  the  execution  rate  of  the  vectorized  code 
on  the  Cray  is  nearly  thirty  times  faster  than  the  rate  on  the  Wright- 
Patterson  AFB  Cyber  170-750.  The  Cray  is  five  times  faster  in  the  scalar 
mode  and  an  additional  factor  of  six  results  from  the  vectorization. 
Because  of  the  large  difference  in  the  time  step  necessary  for  stability 
and  high  frequency  resolution  compared  with  the  fundamental  period  of 
the  oscillation,  the  computation  time  is  extensive.  For  the  full  length 
inlet,  central  processor  time  is  in  excess  of  one  hour  per  buzz  cycle. 
Additionally,  the  core  memory  used  with  the  Cray  is  well  in  excess  of 
that  available  on  the  Cyber.  The  use  of  the  Cray  has  been  an  essential 
element  in  the  successful  completion  of  this  work. 


Section  IV 


DISCUSSION  OF  RESULTS 

4.1  Experimental  Test  Description  and  Results 
Experimental  Apparatus 


The  experimental  work  of  Nagashima,  Obokata  and  Asanuma  (Ref  1) 
was  chosen  as  the  basis  for  numerical  comparison.  At  the  time  the 
computational  effort  was  begun,  this  appeared  to  offer  the  most  complete 
data  set  available.  In  order  to  understand  the  numerical  results,  it 
is  necessary  to  consider  in  some  detail  the  experimental  test. 

An  external  compression  axi-symmetric  inlet  was  tested  at  zero 
and  three  degrees  angle  of  attack  in  a  supersonic  blow-down  wind 
tunnel.  The  test  section  Mach  number  was  fixed  at  a  value  of  two. 

From  the  reported  data,  settling  chamber  total  pressure  was  determined 
to  be  three  atmospheres.  Total  temperature  was  assumed  to  be  the 
standard  reference  temperature  of  15eC.  Reynolds  number  based  on  the 
model  diameter  of  6  cm  was  then  set  at  2.35  x  10^.  Run  time  was 
350  seconds  at  three  atmospheres. 

A  cross  sectional  view  of  the  inlet  model  is  presented  in  Figure 
5.  The  length  of  the  inlet  cowl  was  63.5  cm  with  an  L/D  ratio  of  roughly 
15.88.  The  centerbody  half  angle  was  25°.  Three  different  combinations 
of  spacers  were  used  to  vary  -the  position  of  the  centerbody  with 
respect  to  the  cowl.  These  are  identified  in  Ref  (1)  as  spacers  A,  8 
and  C.  Only  spacer  A  was  considered  for  numerical  comparison  because 
the  majority  of  the  reported  data  was  for  this  case.  The  forebody 
geometry  for  spacer  A  is  shown  in  Figure  6  and  Table  1  tabulates  the  geo¬ 
metric  coordinates  of  the  centerbody  and  cowl.  A  D.C.  motor  was  used  to 
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Fig.  6  Inlet  Model  Characteristics  (from  Ref  1) 
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Outer  Cowl  Wall  Geometry  (mm) 
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Table  1  Inlet  Model  Coordinates 


control  the  position  of  a  throttle  valve  In  order  to  vary  the  exit 
area.  This  area  could  then  be  varied  continuously  within  the  throttle 
ratio  limits. 


TR  -  A  /A  -  0  -*>  2.41  (4.1) 

e  c 

where 

Afi  -  upstream  capture  area  at  cowl  lip 

Ag  -  exit  area  of  holes  located  along  side  walls  of  cowling 
(not  the  minimum  or  choking  area) 

Instrumentation 

The  model  was  Instrumented  with  seven  static  pressure  probes 
of  strain  gauge  type.  Two  were  placed  on  the  cone  surface  (p^,  p,,) 
and  one  at  the  throat  (p^) ,  Four  additional  probes  were  placed  with¬ 
in  the  diffuser,  before  (p^,  p^)  and  after  (p^)  the  diffuser  strut, 
and  in  the  plenum  chamber  (p^) .  Probe  locations  are  shown  in  Figures 
5  and  6.  Flow  visualization  of  the  exterior  flow  field  was  provided 
by  high  speed  schlleren  photography  at  5600  frames/sec.  A  power 
spectral  analyser  was  used  to  provide  Fourier  analysis  data  of  the 
unsteady  pressure  traces.  No  other  flow  field  information  was  avail¬ 
able. 


Experimental  Results 

Spacers  A,  B,  and  C  were  tested  for  a  variety  of  throttle  ratios 
corresponding  to  supercritical,  critical,  and  subcritical  flow  condi¬ 
tions.  Attention  will  be  restricted  to  spacer  A  with  which  numerical 
comparisons  will  be  made.  For  this  configuration,  the  angle  between 
the  centerbody  apex  and  the  cowl  tip  was  31.6°.  This  value  correspond 
ed  to  a  design  Mach  number  of  near  four,  at  which  the  oblique  conical 


shock  impinges  upon  the  cowl  lip. 

At  the  test  condition,  Mach  two,  the  flow  was  reported  to  he 
steady  for  throttle  ratios  (TR)  greater  than  1.14.  Pressures  cor¬ 
responding  to  several  supercritical-critical  conditions  are  shown  in 
Figure  7.  Pressure  rise  within  the  duct  was  reported  by  the  authors 
to  be  smooth,  indicating  the  presence  of  a  shock  train  due  to  shock¬ 
boundary  layer  interaction  rather  than  a  simple  normal  shock.  A 
schlieren  picture  of  the  external  shock  structure  is  available  in 
Figure  8. 

Further  reduction  in  the  throttle  area  Initiated  the  onset  of 
the  buzz  phenomenon  associated  with  the  subcritlcal  regime.  This  is 
characterized  by  violent  pressure  oscillations  within  the  duct  and  large 
amplitude  shock  movement  on  the  centerbody.  Typical  pressure  traces 
are  shown  in  Figure  9  and  the  corresponding  schlieren  pictures  are 
presented  in  Figures  11  and  12. 

Referring  to  Figures  9  and  10,  there  were  apparently  two  distinct 
buzz  regimes.  At  high  values  of  subcritlcal  throttle  ratio,  the  pressure 
traces  exhibited  a  dominant  frequency  near  the  fundamental  mode  upon 
which  much  higher  frequency  harmonics  were  superposed  in  an  erratic 
manner.  The  shock  wave  movement  and  the  pressure  traces  at  the  front 
and  back  ends  of  the  inlet  were  step-like  in  nature.  Interior  to  the 
duct,  the  wave  form  was  more  sinusoidal.  This  organized  pattern  did 
not  hold  near  the  throat  where  the  pressure  trace  became  highly  irregular. 
The  lower  leg  of  the  bow  shock  was  also  observed  to  oscillate  with  a 
low  amplitude  high  frequency  motion  secondary  to  its  primary  motion. 

As  the  throttle  ratio  was  further  decreased,  the  dominant  fre¬ 
quency  suddenly  shifted  to  a  value  three  times  the  fundamental  with 


Fig.  9  Experimental  Pressure  Fluctuation  (from  Ref  1) 


Fig.  10  Amplitude  and  Frequency  Characteristics  of 
Pressure  Oscillations  (from  Ref  1) 


the  6th  and  9th  inodes  also  observed  from  spectral  analysis.  This  held 
true  until  the  throttle  valve  was  completely  closed,  where  tv  Tunda- 
mental  frequency  was  also  observed  at  the  two  downstream  probes. 
Acoustic  theory,  outlined  in  the  next  section,  predicts  only  odd  fre¬ 
quency  modes  for  a  duct  with  an  open  and  a  closed  end.  The  amplitude 
of  the  fluctuations  at  first  decreased  and  then  reached  a  much  higher 
value  for  the  low  subcrltical  throttle  ratios.  The  shock  movement  was 
larger  in  amplitude  and  appeared  to  be  a  continuous  motion. 


4.2  Theoretical  Considerations 


Before  considering  the  numerical  results,  valuable  insight  into 
the  physical  mechanisms  which  govern  the  buzz  cycle  can  be  gained  by  a 
brief  synopsis  of  related  theoretical  work.  Rockwell  and  Nandasher 
(Ref  65)  and  Hankey  and  Shang  (Ref  66)  both  have  outlined  the  neces¬ 
sary  conditions  for  a  self -sustained  oscillation  to  occur: 

(a)  upstream  feedback  of  reflected  disturbances  from  a  down¬ 
stream  impingement  surface  to  an  unstable  shear  layer  in 
proper  phase  and  at  a  frequency  where  amplification  is  possible 

(b)  selective  amplification  of  disturbances  in  the  unstable 
shear  layer  over  a  limited  frequency  range. 


Feedback  Mechanism 

The  dominant  feedback  mechanism  for  the  inlet  is  the  reflection 
of  expansion  and  compression  waves  in  like  sense  from  the  downstream 
throat.  The  time  required  for  a  wave  to  travel  from  the  upstream 
disturbance  to  the  downstream  throat  and  back  is 


L  L 

"  f  (1+M)c  dX  +  f  (1-M) 


(4.2) 


where 


c  -  local  acoustic  speed  at  which  the  disturbance  propagates 


relative  to  the  mean  flow 
M  -  undisturbed  mean  flow  Mach  number  . 

The  buzz  cycle  begins  with  the  forward  expulsion  of  the  bow  shock  on 
the  centerbody  and  a  consequent  expansion  wave  which  propagates  down 
the  duct.  The  reflected  expansion  wave  interacts  with  the  bow  shock 
and  the  shock  retreats  inward.  This  generates  a  compression  wave 
which  similarly  travels  downstream  and  is  reflected  to  expel  the  bow 
shock  and  begin  the  cycle  again.  If  it  is  further  assumed  that  the  flow 
is  uniform  along  the  length  of  the  duct,  L,  then  the  fundamental  period, 
upon  integration,  is 
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P  -  2T  «  0 

P  (l-M^c 


(4.3) 


and  the  fundamental  and  higher  odd  frequency  modes  as  predicted  for  a 
duct  with  an  open  and  a  closed  end  are 


f  n  -  f  (2n-l)  -  ^  (1-M2)  (2n-l)  n  -  1,2,3,... 


(4.4) 


This  very  simple  model  serves  to  explain  several  of  the  features  observed 
in  buzz.  The  step-like  behavior  of  the  shock  and  the  pressure  probes 
at  the  extreme  ends  of  the  inlet  is  a  result  of  the  time  delay  involved 
in  the  propagation  of  forward  and  rearward  traveling  pressure  waves. 
Equation  (4.4)  also  serves  to  explain  the  experimentally  observed  shift 
in  the  fundamental  frequency  with  throttle  ratio.  As  the  throttle 
valve  is  closed,  the  average  duct  Mach  number  must  decrease  since  the 
throat  remains  choked  and  the  upstream  flow  conditions  are  determined 
by  the  area  ratio,  (A^/A^) .  For  example,  at  the  condition  TR-.97, 
the  experimentally  observed  frequency  was  f.  ■  110  H  .  This  requires 
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an  average  duct  Mach  number,  -  .36,  which  is  entirely  reasonable. 
Higher  frequency  modes  may  be  interpreted  as  multiple  waves  within  the 
duct. 


Shear  Layer  Instabilit 


The  central  role  played  by  the  separated  boundary  layer  region  has 
been  identified  from  the  earliest  experimental  investigations.  The 
separation  was  found  to  be  produced  either  by  the  slip  line  resulting 
from  the  intersection  of  the  bow  and  conical  shocks  (Ferri  condition. 

Ref  5)  or  as  9  result  of  separation  on  the  centerbody  induced  by  the 
bow  shock  (Dailey  condition,  Ref  10) .  For  the  present  off-design  condi¬ 
tion  the  Ferri-condition  is  not  applicable  since  the  slip  line  does  not 
enter  the  cowl.  Dailey  further  postulated  that  the  separated  region 
would  instantaneously  block  the  inlet  flow  as  it  became  very  large. 

Separated  flows  possess  a  point  of  inflection  in  the  velocity  pro¬ 
file,  u  -0.  Such  flows  have  been  known  to  be  unstable  since 
ynyn 

Rayleigh  (Ref  67)  first  proved  "the  existence  of  a  point  of  inflection 
as  a  necessary  condition  for  the  occurrence  of  instability,''  (Rayleigh's 
first  theorem).  Spatial  stability  theory  has  been  used  by  a  number 
of  authors  to  predict  the  frequency  at  which  the  shear  layer  will 
amplify  small  disturbances.  Michalke  (Ref  68)  and  later  Roscoe  and 
Hankey  (Ref  69)  considered  the  hyperbolic  tangent  velocity  profile.  In 
a  subsequent  study,  Verma,  Hankey,  and  Scherr  (Ref  70)  studied  profiles 
corresponding  to  the  lower  branch  of  the  Falkner-Skan  equations.  The 
approach  was  to  linearize  the  governing  equations  about  some  known 
parallel  flow  solution  for  small  harmonic  disturbances  which,  in  the 
inviscld  limit,  resulted  in  the  Rayleigh  equation.  The  eigenvalue  problem 
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was  then  solved  to  predict  the  frequency  range  where  disturbance 
amplification  was  possible.  The  frequency  range  at  which  amplifica¬ 
tion  was  a  maximum,  or  "natural  frequency",  then  corresponded  to  the 
most  probable  oscillation  frequency.  For  a  self-sustained  oscillation 
to  develop,  the  feedback  must  occur  in  proper  phase  and  at  a  frequency 
where  amplification  is  possible.  The  actual  frequency  of  oscillation 
is  then  fixed  at  the  nearest  resonant  mode,  equation  (4.4),  where  am¬ 
plification  is  positive. 

Although  stability  theory  is  a  powerful  tool  for  understanding 
shear  layer  instability,  it  is  of  limited  value  as  a  predictive  tool. 
Instability  growth  is  bounded  by  viscous  dissipative  effects  and  a 
limit  cycle  is  reached  which  is  not  predicted  by  linear  theory.  For 
this  reason,  the  theory  can  not  predict  the  amplitude  of  the  oscilla¬ 
tions.  A  further  problem  is  that  the  velocity  profile  of  the  separated 
shear  layer  at  the  time  of  the  disturbance  must  be  known  before  a 
stability  analysis  is  possible.  Normally  this  is  not  known  unless  the 
solution  to  the  problem  itself  is  known.  Also,  complex  flowfields, 
such  as  the  subcritical  inlet,  may  exhibit  a  number  of  unsteady  separat¬ 
ed  regions  for  which  a  time-independent  mean  flow  profile  can  not  be 
specified.  Numerical  solutions  to  the  full  governing  equations  are  re¬ 
quired  to  adequately  describe  the  fluid  dynamic  behavior. 

4.3  Numerical  Results 

Overview  of  Computational  Results 

Numerical  calculations  were  peT  formed  over  the  range  of  conditions 


tested  experimentally.  One  super-critical  (near  critical)  case,  TR  •  1 
was  chosen  because  of  its  Intrinsic  Interest  and  also  to  validate  the 
code  before  attempting  the  unsteady  cases.  Subcritical  flow  was 


investigated  at  both  high  (TR  -  ,97)  and  low  (TR  -  0.0)  values  of  inlet 
mass  flow.  Because  computer  time  is  proportional  to  the  period  of  the 
oscillation  and  hence  to  the  length  of  the  inlet  (Equation  4.3),  a 
shorter  configuration  was  considered  before  the  full  length  inlet  was 
attempted.  The  shorter  inlet  was  obtained  from  the  full  scale  model 
by  removing  the  segment  of  length  corresponding  to  the  constant  center- 
body  radius  section  from  25  to  47.7  cm.  This  then  corresponded  to  an 
inlet  with  an  L/D  ratio  of  10.23. 

Geometry  for  Numerical  Models 

The  inlet  geometry  of  the  full  scale  configuration  was  identical 
with  the  experimental  model  with  several  exceptions  which  were  made  to 
facilitate  the  numerical  calculations  without  significantly  affecting 
the  fluid  dynamics.  The  struts  which  position  the  centerbody  in  the 
inlet  were  not  and  can  not  be  modeled  unless  the  full  three  dimensional 
equations  are  solved.  Furthermore,  the  centerbody  was  extended  to 
merge  with  the  throttle  valve,  maintaining  a  minimum  radius  of  .3  cm. 
This  avoided  the  necessity  of  having  to  provide  additional  axial  re¬ 
solution  of  the  wake  which  resulted  when  the  centerbody  was  terminated. 
It  also  avoided  the  potential  difficulty  of  numerically  resolving  the 
region  very  near  the  coordinate  system  singularity  with  a  highly 
stretched  grid.  The  extended  centerbody  preserved  98%  of  the  crossec- 
tional  area  in  that  region  and  had  a  negligible  effect  on  the  total 
volume  of  the  diffuser. 

The  most  important  model  difference  concerned  the  downstream 
throat.  For  the  conditions  tested,  this  throat  remained  choked  at  all 
times.  However,  as  shown  in  Appendix  E,  the  choking  area,  upon  which 
the  sonic  surface  must  terminate,  is  not  the  exit  area,  Ae,  used  in 


equation  (4.1).  It  is,  rather,  the  area  obtained  by  revolving  about 

the  symmetry  axis  the  normal  to  the  throttle  valve  which  intersects 

the  rear  edge  of  the  cowl.  The  two  areas  are  related  if  the  geometry 

of  the  throttle  valve  and  the  blockage  due  to  the  presence  of  the  struts 

in  the  exit  are  known.  Since  this  latter  information  was  not  given  in 

the  report,  a  procedure  was  developed  to  find  the  actual  choking  area, 

A*,  given  the  throttle  ratio  and  hence  the  exit  area,  A  ,  such  that  the 
*  e 

resulting  pressure  ratios  agreed  with  those  of  Figure  7.  The  procedure 
is  outlined  in  Appendix  E,  This  corrected  throttle  ratio  is  then 
referred  to  as  the  area  ratio,  AR, 

AR  -  A. /A  (4.5) 

*  c 

where 

A±  -  minimum  downstream  exit  area  upon  which  the  sonic  surface 
must  terminate. 

The  choking  condition  is  a  sufficient  constraint  for  the  downstream 
flow,  both  physically  and  numerically.  It  suffices  to  position  the 
shock  within  the  duct,  for  given  inflow  conditions,  such  that  mass  con¬ 
tinuity  is  preserved  for  the  resulting  distribution  of  shock  and  viscous 
total  pressure  losses.  In  an  unsteady  flow,  the  choking  condition  will 
correctly  reflect  rearward  propagating  compression  and  expansion  waves. 

If  a  divergent  section  is  added  downstream,  the  flow  may  be  accelerated 
to  a  supersonic  exit  condition.  The  no-change  condition  (Sec  3.3)  may 
then  be  implemented  as  a  numerical  boundary  condition  which  has  no  effect 
upon  the  flow  upstream  of  the  throat.  This  arrangement  provides  a  natural 
way  to  accomodate  the  large  amplitude  unsteady  fluctuations  which  exit 
through  the  downstream  boundary  in  the  subcritical  regime. 

Upon  incorporating  the  changes  outlined  above,  the  geometries 


corresponding  to  the  full  scale  and  short  inlets  computed  numerically 
are  shown  in  Figures  13  and  22.  In  the  present  approach,  the  throttle 
valve  can  not  be  continously  varied,  as  was  done  experimentally,  but 
must  be  changed  impulsively. 

Near-Critical  Steady  State  Solution 

(TR"1.42,  AR-1.16,  L/D-15.88) 

One  steady  state  condition  was  chosen  for  computation.  Primary 
attention  was  focused  upon  the  proper  treatment  of  boundary  conditions, 
mesh  spacing,  turbulence  modeling,  the  extent  of  the  computational 
domain  and  graphical  display.  In  this  manner,  confidence  in  the 
solution  procedure  was  established  before  attempting  the  unsteady 
calculations. 

The  case  chosen  for  comparison  corresponded  to  a  throttle  ratio, 
TR-1.42,  which  is  equivalent  to  an  area  ratio,  AR-1.162.  The  compu¬ 
tational  grid  used  is  shown  in  Figure  14.  There  were  190  points  in 
the  axial  direction  and  64  points  in  the  radial  direction  of  which  30 
were  within  the  duct.  Axial  resolution  was  concentrated  in  forebody 
and  nozzle  regions  of  the  inlet  where  large  flow  gradients,  including 
the  presence  of  shocks  and  separated  boundary  layers,  were  expected. 
Minimum  axial  step  size  was  0.1  cm  and  the  maximum  was  0.6  cm.  Minimum 
radial  step  size  at  the  wall  boundaries  was  set  by  requiring  a  uni¬ 
form  step  size  across  the  duct  height  at  the  cowl  lip.  The  resulting 
^min  “  0*015  cm  was  then  maintained  along  the  length  of  the  inlet. 

This  minimum  step  size,  although  necessary  for  reasonable  computer 
time  with  an  explicit  code,  resulted  in  somewhat  marginal  viscous 
resolution.  For  example  the  minimum  step  size  was  roughly  twice  the 


Table  1  Static  Pressure  Probe  Locations 
(L/D  -  15.88,  see  Fig.  13) 
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theoretical  boundary  layer  height  on  the  centerbody  at  the  cowl  tip. 
Further  downstream  within  the  inlet,  the  turbulent  non-dimensional 
step  size  was 
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which  is  at  the  edge  of  the  viscous  sublayer  of  the  turbulent  bound¬ 
ary  layer.  Although  skin  friction  and  heat  transfer  rates  could 
not  be  predicted  accurately,  the  overall  flow  behavior,  including 
the  location  and  extent  of  separated  regions,  was  presumed  to  be 
computed  with  sufficient  resolution. 

Figure  15  presents  a  Mach  contour  plot  of  the  entire  inlet 
flowfield.  Expanded  views  of  the  forebody  and  nozzle  regions  are 
shown  in  Figure  16.  The  conical  forebody  may  be  compared  with  the 
Taylor-Maccoll  inviscld  solution  for  the  cone.  At  Mach  two,  the 
25°  cone  generated  a  conical  shock  at  an  angle  of  42.5°.  The  flow 
behind  the  shock  was  then  constant  on  rays  emlnating  from  the  cone  tip. 
The  computed  solution  was  in  close  agreement  with  the  shock  angle, 
surface  pressure  ratio  and  the  cone  surface  Mach  number  exterior  to 
the  boundary  layer 
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A  second  oblique  shock  was  generated  by  the  upper  surface  of  the 
cowl  lip  which  intersected  with  the  bow  shock  at  some  distance  above 
the  cowl  lip.  This  intersection  resulted  yet  another  shock  above 
the  intersection.  A  slip  line,  across  which  magnitudes  of  velocity 
are  discontinuous,  also  originated  from  the  intersection  since 
adjacent  streamlines  passing  through  the  two  different  shock  systems 


Near-Critical  Steady  State  Solution 


Fig.  16  Expanded  Views  of  Forebody  and  Downstream  Throat 
Regions  (TR=1.42,  AR=1.16) 


have  the  same  pressures  and  flow  direction  but  different  values  of 
entropy  and  total  pressure.  Further  downstream,  an  expansion  fan  was 
generated  as  the  flow  was  turned  about  the  convex  corner  on  the  outer 
cowl  surface.  As  a  result  of  the  Interaction  with  the  expansion,  the 
outer  shock  was  then  curved  and  the  flow  downstream  of  the  interaction 
rotational.  The  exterior  flow  fields  of  the  computational  and  ex¬ 
perimental  results  may  be  compared  in  Figures  8  and  16. 

The  flow  interior  to  the  duct  was  found  to  be  strongly  dependent 
upon  the  assumptions  made  concerning  the  nature  of  the  flow,  as  to 
whether  it  was  laminar  or  turbulent,  where  the  transition  region 
occurred,  and  how  the  turbulent  eddy  viscosity  was  modeled.  Unfor¬ 
tunately,  there  was  no  experimental  information  to  aid  in  this  respect 
Neglecting  any  tunnel-produced  freestream  turbulence,  the  flow  was 
assumed  to  be  initially  laminar.  At  some  point  downstream,  a  critical 
Reynolds  number  was  reached  in  which  the  laminar  flow  became  unstable 
and  began  a  transition  to  a  turbulent  state.  A  number  of  external 
effects  may  influence  the  process  (Ref  71:  Chap  XVII).  However,  the 
large  adverse  pressure  gradient  associated  with  the  normal  shock 
structure  is  by  far  the  dominant  effect,  which  may  be  presumed  to 
trigger  the  transition  rather  abruptly  downstream  from  the  shock.  The 
shock  occurred  at  a  Reynolds  number  based  upon  forebody  running  length 
Refl  »  1.5  x  10*\  Since  this  is  roughly  half  the  value  of  Reynolds 
number  at  which  transition  begins  on  a  cone  of  the  same  half  angle  and 
equivalent  flight  conditions  (Ref  72),  the  assumption  of  laminar  flow 
to  this  point  appears  to  be  valid.  At  the  normal  shock,  the  Cebeci- 
Smith  eddy  viscosity  model  of  Section  2.4  was  implemented  on  both  the 
centerbody  and  cowl  as  an  approximate  representation  of  the  Reynolds 
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stress  terms  in  the  mean  flow  equations.  Transition  was  given  by  Equa¬ 
tion  2.38  with  the  starting  location  (s^)  specified  at  the  shock  and 
the  final  location  (s^)  a  distance  of  3Ax  downstream. 

With  the  assumed  turbulence  model,  the  internal  flow  structure 
appeared  as  in  Figures  15  and  16.  Because  the  internal  cowl  lip 
incidence  angle,  $c  «  19.24°,  was  aligned  with  the  flow  angle  at  the 
tip,  no  shock  was  generated  on  the  lower  surface  of  the  cowl.  The 
normal  shock  structure  within  the  inlet  is  perhaps  the  most  inter¬ 
esting  aspect  of  the  solution.  Within  the  divergent  section  downstream 
of  the  lip,  the  flow  accelerated  to  an  invlscld  core  Mach  number  of 
1.70  and  a  one  dimensional  mass  averaged  Mach  number  of  1.60.  This  was 
then  followed  by  a  normal  shock.  Although,  with  the  contour  levels  dis¬ 
played,  this  appears  in  Figure  16  as  a  single  normal  shock,  the 
numerical  solution  displayed  considerably  more  complexity.  The 
dominant  shock  was  followed  by  a  region  of  roughly  three  duct  heights 
at  the  cowl  lip,  [youtgj.~yjinugj.3np*  *n  which  the  flow  oscillated 
within  the  Mach  number  range  1.1  to  .9.  A  small  separated  region 
developed  on  the  centerbody  and  extended  downstream  for  a  distance 
of  nine  duct  heights  before  reattaching.  The  boundary  layer  on  the 
cowl  was  thinner  and  remained  attached.  This  large  deviation  from 
ideal  shock  behavior  can  be  attributed  to  shock-boundary  layer 
interaction.  Generally,  as  the  boundary  layer  becomes  thicker,  the 
single  normal  shock  is  replaced  by  a  series  of  weak  blfuricated  normal 
shocks  in  which  the  normal  component  becomes  progressively  smaller. 

This  shock  structure,  termed  a  "pseduo-shock"  by  Crocco,  is  illustrated 
in  Figure  17  and  has  been  discussed  by  several  authors,  (Ref  73:  Chap 
28,  Ref  74:  Sec  B.5,  Ref  75).  Because  of  the  relatively  coarse  grid 


spacing,  the  detailed  structure  of  the  weak  pseudo-shock  region  was 
smeared  although  it  was  presumed  to  be  correct  in  an  integral  sense. 
The  "pseudo  shock"  was  then  followed  by  a  region  of  subsonic  diffu¬ 
sion  from  the  sonic  condition  with  significant  viscous  dissipative 
losses. 


Fig.  17  Typical  "Pseudo-Shock"  Structure 
As  an  important  consequence  of  the  deviation  from  ideal  shock 
behavior,  the  pressure  rise  does  not  occur  discontinuosly  but  grad¬ 
ually  through  the  shock  train  and  the  subsequent  subsonic  diffusion. 
Static  pressure  on  the  centerbody  and  cowl  surface  is  plotted  in 
Figure  18.  Experimental  results  from  Figure  7  are  also  shown.  The 
agreement  between  the  two  was  quite  good.  One  discrepancy  was  noted 
in  that  the  shock  was  positioned  slightly  upstream  of  probe  p^  while 
experimentally  it  occurred  somewhat  downstream.  This  may  be  attri¬ 
buted  to  small  Inaccuracies  in  the  computed  total  pressure  loss  due  to 
shock  and  viscous  effects  which  position  the  shock,  and  to  some  in¬ 
accuracy  in  the  conversion  of  throttle  ratio  to  an  equivalent  area 
ratio. 

Near  the  exit,  the  flow  reaccelerated  from  a  low  subsonic  Mach 
number  (0.3-0. 4)  to  a  sonic  condition  at  the  throat.  Figure  16. 

Several  mass-averaged  quanitles  of  interest  are  shown  in  Figures  19, 
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20  and  21.  Figure  19  presents  mass-averaged  total  pressure  along 
the  inlet  length.  The  majority  of  the  total  pressure  loss  occurred 
in  the  "pseudo-shock"  and  subsonic  diffusion  regions  with  smaller 
dissipative  and  frictional  losses  through  out  the  remainder  of  the 
duct.  Figure  20  plots  mass-averaged  Mach  number.  It  should  be 
noted  that  the  Mach' number  jump  predicted  for  an  ideal  normal  shock 
required  a  distance  of  7.6  cowl  lip  heights  over  28  axial  grid  sta¬ 
tions  to  develop  due  to  viscous-inviscid  Interaction.  Figure  21,  in 
which  the  ratio  of  mass  flux  to  supercritlcally  captured  mass  flux 
is  plotted,  demonstrates  the  steady  state  nature  of  the  flow  at  this 
throttle  ratio.  The  two  "kinks”  were  due  to  truncation  error  in  the 
vicinity  of  the  shock  and  the  nozzle  throat.  However,  in  keeping  with 
the  use  of  conservative  variables  and  shock  capturing  methods,  con¬ 
tinuity  was  maintained  globally. 


Subcrltical  Regime  (TR  *  .97,  AR  ■  .83,  L/D  ■  10.23) 

In  the  higher  throttle  ratio  subcrltical  regime,  nearly  all  the 
experimental  data  was  for  the  condition  TR  ■  .97.  This  case  was  then 
chosen  for  numerical  comparison.  As  explained  earlier,  a  shorter 
inlet.  Figure  22,  was  computed  before  considering  the  full  scale 
model.  This  approach  provided  a  modest  saving  in  grid  points  and  a 
substantial  reduction  in  the  fundamental  period  of  the  oscillation, 
Equation  (4.3)  and,  consequently,  computer  time.  In  experimental  studies 
in  which  the  effect  of  L/D  was  studied  (Refs  9,  10),  the  onset  of  the 


instability  was  unchanged  and  the  shock  wave  movement  and  pressure 


trace  history  were  similar  but  occurred  at  a  corresponding  higher 
frequency  as  L/D  was  decreased. 


While  the  numerical  simulation  was  not  successful  in  capturing 
the  naturally  occurring  self-sustained  oscillations  reported  experi¬ 
mentally,  it  did  provide  insight  into  two  important  aspects  of  the 
phenomenon.  First,  the  calculation  was  found  to  be  extremely  sensitive 
to  turbulence  modeling  for  which  no  experimental  guidance  was  available. 
Although  the  adequacy  of  a  simple  algebraic  eddy  viscosity  model  is 
also  questionable,  the  stability  of  the  flow  was  found  to  be  depend¬ 
ent  upon  the  location  and  extent  of  the  transitional  region.  Second, 
numerical  evidence  was  found  to  support  the  contention  that  the 
separated  boundary  layer  instability  and  the  resulting  flow  block¬ 
age  is  the  mechanism  by  which  the  buzz  cycle  is  triggered. 

A  grid  of  dimensions  154  x  64  was  used  in  this  case.  It  was 
similar  to  Figure  14  but  the  external  domain  was  larger  to  contain 
the  anticipated  shock  movement.  As  in  the  steady  state  case,  the 
flow  was  assumed  to  be  turbulent  downstream  from  the  normal  or  bow 
shock  on  both  the  centerbody  and  cowl.  A  starting  solution  was  ob¬ 
tained  from  the  previous  steady  state  solution  with  the  downstream 
throat  area  Impulsively  decreased.  This  generated  a  compression  wave 
which  traveled  forward  to  force  the  normal  shock  out  onto  the  centerbody. 
The  solution  was  continued  and  converged  to  a  subcrltical  steady 
state  in  which  the  duct  pressure  corresponded  to  the  experimentally  de¬ 
termined  maximum  value  before  buzz  onset.  Mass  flux  through  the  inlet 
was  reduced  to  85%  of  that  captured  in  the  supercritical  regime.  Unfor¬ 
tunately,  the  experiment  indicated  that  the  flow  was  unstable  and  buzz 
present  at  this  condition.  The  computed  turbulent  boundary  layer  pro¬ 
file  exhibited  only  a  small  separated  region  on  the  centerbody.  This 
is  normally  what  is  expected,  since  the  turbulent  boundary  layer  is 
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less  prone  to  separate  than  Its  laminar  counterpart  in  the  presence 
of  adverse  pressure  gradients. 

A  key  hypothesis  presumed  a  large  region  of  separated  flow  as 
a  necessary  condition  for  buzz,  (Sec  4.2).  Consequently,  the  transi¬ 
tion  region  was  moved  downstream  beginning  at  6  cm  (Re  -2.3x10^)  and 
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ending  at  12  cm  (Re  -4.7x10  ).  Although  these  locations  were  rather 

®f 

arbitrary,  the  intent  was  to  allow  a  large  region  of  separated  flow  to 
develop  just  downstream  of  the  cowl  lip  where  blockage  might  then  occur. 
The  solution  was  restarted  from  the  previous  solution  near  the  time  at 
which  the  normal  shock  had  been  forced  from  the  inlet. 

As  expected,  the  flow  was  no  longer  steady.  Figure  23  presents  a 
sequence  of  Mach  contour  plots  which  trace  the  flow  development  in  re¬ 
sponse  to  the  shift  in  transition  location.  A  very  large  region  of 
unsteady  separated  flow  developed  on  the  centerbody  downstream  from  the 
shock  and  extended  to  the  arbitrarily  prescribed  location  for  fully 
turbulent  flow  before  reattaching.  Separation  appears  in  the  Mach  con¬ 
tour  plots  of  Figure  23  as  a  series  of  irregular  and  closely-spaced  lines 
of  constant  Mach  number.  In  results  to  be  presented  later  (Figures  28, 

29)  this  can  be  verified  by  direct  comparison  of  Mach  contour  and  velocity 
plots .  Separation  also  occurred  on  the  cowl  and  the  core  flow  followed 
a  sinuous  path  through  this  region.  Figure  24  plots  the  instantaneous 
mass  flux  through  the  inlet  at  times  corresponding  to  Figure  23.  Taken 
together  it  is  possible  to  at  least  partially  confirm  a  hypothesis  made 
by  Dailey  (Ref  10).  Referring  to  Figure  24,  the  large  separated  region 
was  found  to  reduce,  or  block,  the  mass  flow  through  the  inlet  to  a 
value  less  than  half  that  captured  when  the  separated  region  we*  not 
present.  In  direct  response  to  this  blockage,  the  shock  was  forced  to 
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its  most  forward  position  at  frame  t  -  .494x10  sec  before  retreating  to 

_2 

its  most  rearward  position  at  frame  t  ■  .805x10  sec.  It  thus  appears 
that  the  origin  of  the  instability  was  within  the  unsteady  separated  bound¬ 
ary  layer  region  just  inside  the  duct.  The  bow  shock  remained  in  a  stable 
subcritical  position  until  forced  forward  by  the  momentary  downstream 
blockage  associated  with  the  unsteady  motion  of  the  separated  boundary  layer. 

If  the  departure  from  equilibrium  were  of  sufficient  magnitude  the 
inlet  would  then  be  expected  to  enter  the  buzz  cycle.  However,  the 
globally  organized  buzz  motion  did  not  occur.  As  can  be  seen  from  Fig¬ 
ure  23,  the  Instability  remained  relatively  weak  and  localized  near  the 
laminar  cowl  lip  region  and,  in  particular,  did  not  propagate  to  the  aft 
end  of  the  duct.  Although  the  duct  flow  is  almost  certainly  turbulent, 
the  turbulence  model  appeared  to  artificially  damp  the  instability.  Liou 
and  Coakley  (Ref  38)  also  noted  a  substantial  sensitivity  in  the  degree 
separation  and  unsteadiness  to  the  turbulence  model. 

An  additional  difficulty  is  evident  in  Figure  11  in  which  the  ex¬ 
perimental  shock  wave  structure  and  the  flow  field  are  seen  to  be  asym¬ 
metric.  This  was  noted  in  Reference  1  and  attributed  to  slight  center- 
body  ecentrlcity.  Nevertheless,  an  axi-symmetric  solution  was  computed 
since  this  was  not  thought  to  be  a  dominant  mechanism.  It  is  clear  that  the 
asymmetric  shock  structure  would  produce  a  larger  region  of  separated 
flow  on  the  leeward  side  of  the  inlet  and  in  turn  affect  the  downstream 
blockage.  Given  the  sensitivity  of  the  flow  field  simulation  to  the  ade¬ 
quacy  of  the  turbulence  modeling  and  the  complete  lack  of  flow  field 
Information,  the  exact  nature  of  the  instability  process  is  likely  to  re¬ 
main  unknown  until  a  more  detailed  experimental  study  is  undertaken. 

Because  buzz  was  not  found  at  this  condition,  the  full  length  inlet  was 
not  computed. 
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Fig.  23  (cont.) 
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Fig.  24  Ratio  of  Instantaneous  Mass  Flux  to  Supercritically 
Captured  Mass  Flux  (TR*.97,  AR=.83,  L/D«10.23) 
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One  flow  condition  at  very  low  throttle  ratio  was  chosen  for  num¬ 
erical  investigation.  Referring  to  Figure  10,  the  oscillations  ob¬ 
served  at  this  condition  were  considerably  larger  than  at  the  higher 
subcritical  throttle  ratios. 

A  decision  was  made  to  compute  the  case  in  which  the  throttle  valve 
was  completely  closed  for  several  reasons.  From  Figure  10,  it  can  be 
seen  that  the  amplitude  of  the  oscillations  was  similar  to  that  obtained 
when  the  throttle  valve  was  nearly  closed  (TR  «  .35,  .55).  Experimentally, 
at  low  throttle  ratios,  the  frequency  was  observed  to  jump  to  the  third 
and  higher  modes.  This  also  held  true  at  TR  *  0.0  except  that  the  funda¬ 
mental  mode  was  again  evident  at  the  two  downstream  orobes.  The  primary 
reason  for  choosing  this  case,  however,  was  economy.  In  order  to  compute 
the  flow  through  the  very  small  nozzle  throat  at  low  throttle  ratios,  the 
minimum  step  size  across  the  throat  would  be  further  reduced.  Because 
the  code  was  explicit,  the  allowable  time  step  size  would  then  be  set  by 
the  throat  at  a  much  smaller  value  than  was  necessary  for  the  remainder 
of  the  inlet.  The  smaller  time  step  would  have  resulted  in  unacceptable 
increases  in  computer  time.  As  an  alternative,  the  nozzle  configuration 
used  in  the  near-critical  solution  (TR  *  1.42)  was  also  used  in  this  case. 

A  solid  wall  boundary  condition  was  imposed  at  the  throat  station  and 
the  computational  domain  terminated  at  that  point. 

Self-sustained  oscillations  were  found  at  both  L/D  ratios  (15.88, 
10.23)  for  this  throttle  condition.  The  discussion  will  concentrate  on 
the  full  length  inlet  which  corresponds  with  the  experimental  test. 

A  brief  discussion  of  the  shorter  inlet  will  then  be  given. 


The  downstream  grid  for  the  full  length  inlet  was  identical  to 
that  of  Figure  14  with  the  divergent  section  aft  of  the  sonic  throat 
omitted.  The  upstream  region  was  enlarged  to  contain  the  shock  struc¬ 
ture  during  the  buzz  cycle  and  is  shown  in  Figure  25.  Total  grid 
dimensions  were  then  202  x  64.  Initial  conditions  were  fixed  assuming 
stagnation  conditions  within  the  inlet: 

u  ■  v  *>  0 

p  =  maximum  pressure  at  buzz  onset  from  experimental  results 

T  *  adiabatic  wall  temperature 

The  external  flow  field  was  initialized  from  previous  solutions  in 
which  the  shock  had  been  forced  into  a  subcritical  position  on  the 
centerbody. 

Because  the  turbulence  model  had  previously  been  found  to  damp 
the  occurrence  of  the  Instability  it  was  not  applied  in  the  present 
case.  It  was  expected  that  the  unsteady  flow  would  consist  of  a  large 
scale  organized  motion  of  discrete  low  frequency  peaks  and  a  broad  band 
of  higher  turbulent  frequencies.  The  numerical  code  is  capable  of 
directly  calculating  a  finite  number  of  low  frequency  components  up  to 
the  shortest  wave  length  which  can  be  resolved  by  the  grid.  The  tur¬ 
bulence  model  should  then  only  simulate  the  smaller  scale  turbulent 
structure  of  the  unsteady  flow.  For  example,  the  nth  Fourier  component 
of  a  longitudinal  wave  propagating  with  acoustic  speed  may  be  expressed 

88  IK  (x-ct)  I(k  x-o)  t) 

Aen  ■  A  e  n  n  (4.6) 

n  n 

where 

I  -  rCT 

c  “  acoustic  velocity 


k  -  t - wave  number  of  nth  mode 

n  X 

n 

X  -  —  -  wave  length  of  nth  mode 
n  n  6 


The  smallest  wave  length  which  is  spatially  resolved  by  the  mesh  is 
X  *  2Ax,  where  Ax  is  the  largest  step  size  on  the  non-uniform  mesh. 
Thus,  kfflax  *  it/ Ax  and  the  corresponding  frequency  is  established  from 


to  k  c 
f  _n_  n 

n  2ir  2it 

For  the  present  case  (Ax)  ■  0.6  cm  and  f  -  25  KHz. 
r  max  max 


(4.7) 


On  the  other  hand.  Chapman  (Ref  39)  has  quoted  a  Strouhal  number  for  the 
mean  frequency  of  the  turbulent  eddies  for  flow  over  a  flat  plate  as 


st6-r  *0-2 


(4.8) 


If  the  shear  layer  thickness,  6,  is  taken  as  that  of  the  separated 
boundary  layer  thickness  near  the  cowl  lip,  and  the  above  formula 
assumed  valid  then  f  -  30KHZ.  Thus  the  lowest  Fourier  components 

t  *6* 

which  provide  the  principle  transport  of  momentum  and  energy  are 
resolved  while  the  higher  Fourier  components  which  contain  the  dis¬ 
sipative  processes  are  not  resolved  and  should  properly  be  modeled. 

This  implies  that  the  period  over  which  the  fluctuating  quantities  are 
averaged  should  be  grid  and  problem  dependent  (Ref  46) .  Existing  tur¬ 
bulence  models,  derived  from  stationary  data  bases,  may  overpredict 
the  appropriate  eddy  viscosity.  If  the  turbulence  model  is  omitted,  the 
transport  processes  are  resolved  while  the  dissipative  processes  are  not. 


Shang  (Ref  36)  has  recently  shown  that  solutions  to  the  laminar  Navier- 
Stokes  equations  can  reproduce  some  "large-scale"  turbulent  statistical 
properties  for  an  unsteady  oscillatory  flow  in  the  transitional  Reynolds 
number  range.  Although  the  fundamental  frequency  of  the  buzz  motion  is 
on  the  order  of  100  HZ,  the  small  scale  motion  of  the  separated  boundary 
layer  occurs  at  a  much  higher  frequency  This  motion  should  not  be 


supressed  by  the  turbulence  model  in  the  numerical  simulation  of  the 
Instability. 

Under  these  assumptions,  numerical  solutions  were  computed  through 
three  complete  buzz  cycles.  The  instability  was  found  to  develop  im¬ 
mediately  as  a  consequence  of  the  non-equilibrium  state  of  the  initial 
conditions.  Figure  26  presents  pressure-time  traces  at  locations 
corresponding  to  the  seven  experimental  pressure  probes.  Additional 
locations  are  shown  in  Figure  27  for  five  forebody  stations  and 
thirteen  equally  spaced  locations  within  the  duct  on  the  cowl  surface. 

All  probe  locations  are  shown  in  Figure  13.  As  can  be  seen  from  the 
pressure  traces,  the  instability  was  first  evident  near  the  cowl  lip 
region  with  rapid  shock  expulsion  on  the  centerbody  well  before  the 
disturbance  was  felt  downstream.  The  oscillation  was  found  to  quickly 
reach  a  bounded  state  in  which  the  amplitude  was  constant.  The 
oscillatory  motion  may  be  described  as  quasi-periodic  since  the  pressure¬ 
time  traces  do  display  a  definite  period  but  are  not  repeatable  due 
to  the  irregular  superposition  of  higher  frequency  components.  The 
dominant  step-like  behavior  of  the  wave  form  is  obvious  at  all  locations 
except  in  the  region  near  and  immediately  downstream  from  the  cowl  lip. 
The  frequency  of  this  dominant  mode  can  be  accurately  established  by 
evaluating  the  period  at  which  the  shock  just  crosses  any  of  the 
forebody  locations  on  successive  cycles.  The  calculated  frequency 
was  found  to  be  128  HZ.  This  agrees  closely  with  the  simple  theoretical 
prediction.  Equation  (4. A),  of  f^  -  127  HZ.  Additional  higher  frequency 
modes  were  introduced  by  the  unstable  separated  boundary  layers  near 
the  cowl  lip  region  where  no  recognizable  pattern  was  evident.  Since 
only  three  cycles  were  computed,  there  was  insufficient  data  to  perform 
a  Fourier  analysis  of  the  wave  form. 
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A  sequence  of  41  Mach  contour  plots,  defining  the  third  buzz  cycle 


is  shown  in  Figure  28.  The  bow  shock  was  forced  to  the  tip  of  the 


centerbody  as  a  result  of  interaction  with  the  reflected  compression 


wave.  In  the  explusion  phase,  a  region  of  reverse  flow  was  found  to 


extend  between  the  base  of  the  shock  and  the  cowl  lip.  The  strong 


shear  layer  dividing  the  two  flow  regions  is  clearly  visible  in  frames 


t  ■  .159  x  10  ^  -*•  t  *  .171  x  10"^  sec.  In  frame  t  «  .173  x  10  ^  sec, 


the  shock  has  been  pushed  to  the  tip,  at  which  point  the  shear  layer 


ruptured  and  flow  was  then  expelled  radially  outward  from  the  inlet. 


The  shock  then  remained  near  the  centerbody  tip  and  flow  expulsion 


continued  until  frame  t  •  .197  x  10  sec.  At  this  point  the  inlet 


started  to  ingest  mass  again  and  the  shock  began  a  rearward  retreat 


towards  the  cowl.  The  shock  then  assumed  a  position  just  forward  of 


the  cowl  lip  where  it  remained  throughout  the  remaining  half  of  the 


buzz  cycle.  During  this  part  of  the  cycle,  the  presence  of  several 


separation  cells  was  clearly  visible.  The  first  extended  from  the 


base  of  the  shock  and  successive  cells  alternated  between  the  center- 


body  and  the  cowl  out  to  about  10  to  12  cm.  The  core  flow  followed 


a  sinuous  path  through  this  region.  The  cell  length  was  found  to  be 


roughly  3  cm.  There  is  some  evidence  that  the  frequency  (~  2800  HZ) 


associated  with  this  length  scale  through  Equation  4.4  may  correspond 


to  some  of  the  more  random  small  amplitude  oscillations  evident  in 


Figures  22  and  23.  A  16  mm  color  graphics  movie  was  also  made  from 


this  same  contour  sequence.  In  the  movie,  the  lower  leg  of  the  bow 


shock  was  observed  to  undergo  a  small  amplitude  oscillatory  movement  in 


addition  to  the  basic  buzz  cycle  movement.  This  motion  was  apparently 


in  response  to  the  transient  flow  blockage  caused  by  the  unsteady  motion 


of  the  separated  boundary  layers  downstream. 
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Figure  29  consists  of  a  sequence  of  21  velocity  plots  again  cov¬ 
ering  the  third  buzz  cycle.  Many  of  the  previously  described  features 
are  also  evident  here  as  well.  In  particular,  the  reverse  flow  region 
which  forces  the  shock  forward,  the  radial  outflow  of  the  expelled 
mass,  and  the  separation  cells  within  the  inlet  are  obvious.  As  a 
result  of  the  mass  expulsion,  a  vortex  structure  situated  on  the  cowl 
lip  can  also  be  seen  in  frames  t  ■  .167  x  10  ^  .179  x  10  *  sec. 

Figure  30  provides  a  time  history  of  mass  flux  through  the  inlet 
(cowl  lip  to  nozzle  throat)  for  the  computed  three  cycles.  The  propa¬ 
gation  of  disturbances  through  the  inlet  can  be  tracked  by  following 
successive  time  frames  although  this  becomes  more  difficult  as  addition¬ 
al  disturbances  are  introduced.  One  suprising  result  was  the  magnitude 
of  the  mass  flux  through  the  inlet.  It  was  found  that  the  maximum  mass 
flux  in  both  the  expulsion  and  ingestion  phases  was  consistently  near 
the  same  magnitude  as  the  supercritically  captured  mass  flow.  The  sep¬ 
arated  boundary  layer,  depending  upon  the  extent  of  separation,  appeared 
to  control  the  mass  flux  ingested  by  the  inlet.  Since  the  downstream 
throat  was  closed,  once  the  oscillation  reached  a  bounded  limit  cycle, 
the  expelled  mass  flux  must  necessarily  be  equal  to  that  ingested.  Given 
the  bounded  nature  of  the  pressure  oscillations,  this  apparently  held 
true  although  the  first  cycle  shock  expulsion  was  somewhat  weaker  and 
the  second  somewhat  stronger  than  that  shown  for  the  third  cycle. 

The  only  experimental  information  available  for  this  case  appears 
in  Figure  10.  The  dominant  frequency  was  apparently  three  times  the 
fundamental  mode  although  the  fundamental  mode  was  observed  at  the  two 
downstream  probes  (p6,  p7) ,  The  numerical  solution  appeared  to  be  a 
hybrid  of  the  two  buzz  regimes  shown  in  Figure  9.  The  fundamental  mode 
was  dominant  and  there  was  considerable  similarity  in  both  the  waveforms 
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and  the  phase  relationships  with  the  TR  -  .97  case.  Figures  9  and  26 
are  replotted  in  Figure  31  for  side-by-side  comparison.  The  ampli¬ 
tude  of  the  oscillations  from  the  numerical  solutions  more  closely 
resemble  the  experimental  result,  TR  ■  G.O.  The  oscillation  amplitude, 
as  obtained  from  Figure  26,  is  plotted  in  Figure  32  and  compared  with 
experimental  values  from  Figure  10.  Because  the  computed  wave  form 
is  different,  the  point  by  point  correspondence  with  Figure  10  is 
not  particularly  good,  although  the  maximum  and  minimum  bands  are  in 
agreement.  Figure  33  depicts  the  time  averaged  values  of  the  pressure 
traces  given  in  Figures  26  and  27  with  the  2AP  band  superposed  on 


the  average  valve.  Numerically,  the  oscillation  amplitude  steadily 
increased  along  the  duct  length  while  the  experiment  displayed  an 
antinode  at  probe  4  as  well  as  probe  7. 

The  sudden  frequency  shift  as  throttle  ratio  is  decreased  at  low 
mass  flows  has  also  been  observed  in  other  experimental  works.  Dailey, 

(Ref  10),  reported  that  the  higher  frequency  does  not  scale  with  L/D 
but  remains  relatively  constant  as  length  is  changed,  suggesting  a  local 
instability  mechanism.  Trlmpi  (Ref  9)  found  the  high  frequency  motion 
to  intermittently  break  down  and  reform.  Nagashima,  (Ref  1),  was  unable 
to  capture  the  higher  frequency  motion  when  the  centerbody  was  removed. 
Although  several  hypotheses,  such  as  vortex  shedding,  (Ref  10),  or  an 
edge  tone  phenomenon,  (Ref  13),  have  been  advanced  to  explain  the  fre¬ 
quency  shift,  this  mechanism  is  still  not  understood  today.  It  was  specu¬ 
lated  that  the  higher  frequency  might  be  triggered  numerically  by  an 
Initial  condition  consisting  of  several  spatial  waves  along  the  inlet  length 
but  this  has  not  been  tested.  Given  the  complexity  of  the  physical  mech¬ 
anisms  Involved,  numerically  capturing  the  fundamental  mode  is  thought  to 
be  an  encouraging  first  step. 
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Fig.  26  Computed  Pressure  Fluctuations  for  Experimental  Probe 
Locations  (see  Fig.  13) ,  (TR«0.0,  L/D=15.88) 
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Fig.  30  Ratio  of  Instantaneous  Mass  Flux  to  Supercritically 
Captured  Mass  Flux  (TR-0.0,  L/D*15.88) 
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The  buzz  cycle  obtained  for  this  length  was  entirely  similar  to  that 


obtained  with  the  full  length  Inlet.  Figure  34  gives  pressure-time 
traces  at  the  nine  probe  locations  shown  In  Figure  22.  As  expected,  the 
fundamental  frequency  increased  to  188  HZ  which  is  in  good  agreement  with 
the  predicted  frequency.  Equation  (4.4),  f^  -  192  HZ.  The  small  amplitude 
higher  frequency  motion  appears  to  be  similar  to  that  found  in  the  full 
length  inlet  and  does  not  scale  with  the  length  of  the  inlet. 
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L  ■  44.35  cm 

TABLE  3  Static  Pressure  Probe  Locations 
(L/D  -  10.23,  see  Fig.  22) 
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Fig.  34  Computed  Pressure  Fluctuations  for  9  Probe 
Locations  (see  Fig.  22),  (TR-O.O,  L/D-10.23) 
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Section  V 


CONCLUSIONS,  ACCOMPLISHMENTS,  AND  RECOMMENDATIONS 

Conclusions 

A  numerical  method  has  been  used  to  solve  the  unsteady,  compres¬ 
sible  Navier-Stokes  equations  for  the  flow  field  about  an  external 
compression  axi-symmetric  inlet  with  a  length  to  diameter  ratio,  L/D  = 
15.88,  at  Mach  2.0,  and  a  Reynolds  number  based  upon  diameter  Re^  ■ 
2.36x10**.  Solutions  were  obtained  in  the  near-critical  and  subcritical 
flow  regimes.  The  near-critical  solution  reached  a  stable  steady  state 
while  the  subcritical  solutions  were  found  to  reach  an  unsteady  bounded 
oscillatory  state  known  as  inlet  buzz. 

The  steady  state  near-critical  solution  agreed  well  with  all  avail¬ 
able  experimental  data.  Investigations  at  the  high  subcritical  throttle 
ratio,  TR  *  .97,  demonstrated  a  large  sensitivity  to  turbulence  modeling 
with  the  model  apparently  tending  to  damp  the  instability.  Evidence  was 
also  found  to  support  the  hypothesis  that  the  separated  boundary  layer 
and  its  flow  blocking  effect  plays  a  key  role  in  the  instability  mech¬ 
anism.  At  low  subcritical  throttle  ratios,  TR  ■  0.0,  buzz  was  computed 
for  two  different  L/D  ratios.  In  both  cases,  the  dominant  frequency 
corresponded  to  the  fundamental  mode  predicted  by  a  simple  wave  propaga¬ 
tion  model,  and  the  experimentally  observed  shift  to  higher  frequency 
modes  was  not  found  numerically.  The  computed  solutions  displayed  the 
large  amplitude  pressure  oscillations  and  traveling  shock  waves  charac¬ 
teristic  of  inlet  buzz.  In  addition,  large  regions  of  separated  flow 
exhibiting  a  high  frequency  instability  together  with  the  feedback 
mechanism  of  reflected  acoustic  waves  were  clearly  shown  in  the  numerical 


solutions . 


Accomplishments 


The  significant  accomplishments  resulting  from  the  present  work 

Include  the  following: 

1.  For  the  first  time,  numerical  solutions  of  the  fundamental  equations 
have  been  used  to  capture  the  naturally  occurring  self-sustained 
oscillations  of  an  inlet  operating  in  the  subcritical  flow  regime. 

2.  The  computed  solutions  provide  further  insight  into  buzz  phenomenon. 

In  particular,  the  critical  importance  of  the  separated  shear 
layer  as  an  amplifier  of  small  disturbances  coupled  with  the 
closed- loop  feedback  of  reflected  disturbances  has  been  shown.  Inlet 
buzz  is  thus  generically  related  to  other  flows  with  self -sustained 
oscillations  and  can  be  analysed  in  that  theoretical  framework. 

3.  The  computed  near-critical  steady  state  solution  was  found  to  be  in 
good  agreement  with  the  experimental  result.  In  particular,  the  gross 
behavior  of  the  "pseudo-shock"  and  subsequent  subsonic  diffusion 
regions  resulting  from  the  boundary  layer  interaction  with  the  terminal 
normal  shock  was  predicted.  To  the  author's  knowledge,  this  is  the 
first  numerical  solution  in  which  this  phenomenon  has  been  observed. 

4.  The  present  work  represents  the  first  use  of  a  choking  throat  in  a 
diffuser  calculation  to  position  the  shock  train  within  the  inlet. 

This  appears  to  be  an  attractive  alternative  to  other  methods  such 
as  specifying  pressure  at  a  downstream  subsonic  boundary.  For  un¬ 
steady  oscillatory  flows  this  approach  allows  the  passage  of  large 
amplitude  disturbances  through  the  downstream  boundary  with  no 
numerical  boundary  condition  influence  on  the  solution. 

In  the  course  of  obtaining  the  above  numerical  solutions; 

5.  An  existing  two  dimensional  code  has  been  converted  to  axi-symmetric 


coordinates  in  strong  conservation  law  form.  The  resulting 
code  was  then  vectorized  for  efficient  use  on  the  CRAY-1S  at  the 
A.F.  Weapons  Laboratory  where  solutions  were  obtained  by  remote 
access. 

6.  .  A  suitable  mesh  was  generated  and  the  algorithm  appropriately 

modified  to  solve  for  both  the  internal  and  external  flow  about  an 
axi-symmetric  inlet. 

Recommendations 

In  the  course  the  present  effort,  several  problem  areas  were  dis¬ 
covered  which  should  be  addressed  in  future  efforts.  The  need  for  a 
more  economical  means  of  obtaining  numerical  solutions  quickly  became 
apparent.  The  solutions  generated  by  the  present  explicit  code  were  only 
possible  on  a  machine  such  as  the  CRAY.  Implicit  codes  with  unconditional 
stability  offer  the  very  attractive  possibility  of  much  larger  allowable 
time  steps.  Balanced  against  this,  the  time  step  size  necessary  for 
accurate  temporal  resolution  is  still  restrictive  and  the  execution  time 
per  Iteration  is  larger.  Also,  the  robustness  for  large  amplitude  os¬ 
cillatory  flows  as  well  as  the  efficiency  on  a  vector  processor  for  such 
a  code  would  need  to  be  verified. 

A  second  need  is  for  a  clearer  understanding  of  the  role  of  turbu¬ 
lence  and  turbulence  modeling  in  the  solution  process.  The  transition 
region  should  be  known  from  experimental  results.  The  simple  algebraic 
eddy  viscosity  model  has  known  deficiencies  for  separated  flows  with 
large  adverse  pressure  gradients.  Turbulence  models  for  such  flows  are 
an  area  of  current  research.  At  present,  it  is  not  clear  that  the  more 
complex  non-equilibrium  turbulence  models  are  superior  for  this  type 
of  flow.  The  eddy  viscosity  model  which,  in  effect,  lowers  the  effective 


Reynolds  number  of  the  flow  by  several  orders  of  magnitude  tended 
to  numerically  damp  the  flow  instability.  Since  the  grid  was  capable 
of  resolving  directly  the  lower  turbulent  eddy  frequencies,  the 
turbulence  model  should  properly  account  for  only  the  higher  fre¬ 
quencies.  It  appears  that  existing  turbulence  models,  derived  from 
stationary  data  bases,  may  overpredict  the  eddy  viscosity  for  unsteady 
mean  flows  with  frequencies  lower,  but  not  orders  of  magnitude  lower, 
than  the  typical  turbulent  eddy  frequencies.  Appropriate  modification 
of  existing  turbulence  models  for  this  class  of  flow  should  be  Investigated 
A  systematic  approach  to  this  problem,  large  eddy  simulation,  is  still 
in  its  infancy. 

A  final  need  is  for  more  detailed  experimental  work.  Good  flow 
visualization  and  detailed  measurements  are  invaluable  in  verifying 
numerical  results.  This  is  particularly  true  for  items  such  as  turbu¬ 
lence  which  must  be  modeled  yet  appear  to  have  a  large  bearing  upon  the 
solution.  In  addition  to  errors  resulting  from  turbulence  modeling,  the 
computed  solutions  are  subject  to  numerical  error.  In  principle,  for 
steady  state  problems,  truncation  error  can  be  removed  by  successive  grid 
refinement  and  techniques  such  as  Richardson  extrapolation.  However,  at 
present,  the  rigorous  application  of  such  procedures  is  prohibitive  due 
to  cost  for  many  realistic  problems.  For  unsteady  flows,  spatial  trunca¬ 
tion  error  varies  as  a  function  of  time  and  grid  refinement  studies 
would  have  to  be  compared  at  corresponding  time  increments.  Unsteady 
problems  are  further  complicated  by  numerical  dispersion  and  phase  errors. 
Experimental  measurements  provide  an  essential  check  on  numerical  solutions 
for  complex  flows  in  which  the  fluid  mechanics  are  not  well  understood. 

In  fact,  numerical  and  experimental  Investigations  are  complimentary  and 


their  joint  use  should  produce  a  synergistic  effect  which  offers  the 
best  hope  of  fully  understanding  the  very  complex  fluid  dynamic  pro¬ 
cesses  underlying  inlet  buzz. 
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Appendix  A 


DIFFERENCE  OPERATORS  USED  IN 
MACCORMACK'S  ALGORITHM  (PHYSICAL  VARIABLES) 

First  Derivatives 

Vx  fi,J  ■  «i,J  -  £i-l.J>£4*  ■  +(-  I<8xxf)i,jta  +  • 


4,  fi,i  ‘  <£1+1,J  -  ■  (3*f)l,J  +{2  + 


sx  fi,j  -  i(,x +  V*i.j  ■  -  £i-i,j)/2ta 

■  <9x£)l,j  +  ,l  <W>i,Jto2  +  -> 

Second  Derivatives 

6xx  fi,j  *  Vxfi,j  *  Vxfi,j  *  (fi+l,J  "  2fi,j  +  fi-l,j)/Ax2 

■  <3xxf)itl  +  (l2  ‘’xxxx^i.J  a*2+  -l 
Sxy  £i,j  ■  T  <Vx  +  V  Sy  fi,j  ■  SxSy  £i,3 

(fi+l,j+l  *  fi+l,j-l  *  fi-l,j-l  “  fi-l,j+l)/4AxAy 

-  O  f)4  x  +{tO  f)4  4Ax2+ f).  . Ay2+ . . . } 
xy  i, j  6'  xxxy  'i,j  6  xyyy  'i,j 

Similar  operators  may  be  defined  for  derivatives  with  respect  to  y,  as 
well  as  derivatives  with  respect  to  the  transformed  variables*  £,n. 


ft  &1- 


Appendix  B 

TRUNCATION  ERROR  FOR  MACCORMACK’S  ALGORITHM 
Because  of  the  asymmetric  differencing  in  time  and  space  it  is  not 
immediately  obvious  what  the  order  of  the  truncation  error  terms  are  in 
MacCormack's  method.  It  is  useful  to  demonstrate  the  second  order  acc¬ 
uracy  of  the  algorithm  rigorously.  In  what  follows,  it  is  not  necessary 
to  consider  the  individual  terms  of  the  F  vector.  In  particular,  dif¬ 
ferencing  of  the  viscous  derivative  terms  can  be  shown  to  be  of  equiva¬ 
lent  accuracy  as  the  lnviscid  terms  and  need  not  be  considered  in  detail 
Also,  because  of  the  use  of  splitting,  it  is  only  necessary  to  consider 
one  spatial  dimension. 


sa  .. 


Given  the  non-linear  system  of  equations 

If  <°>  +  &  <r)  ’ 0 


(B.l) 


where  U(x,t),  F(U)  are  n-dlmensional  vectors  (for  the  present  case,  n-A) 


A  difference  equation  approximation  may  be  written  as 

Ui+1  *  u"  ~  o 

u^1  -  \  [u"  +  U®+1  -  o  (ijJJ  -  fJ+1)] 

where 


(B.2) 


(B.3) 


Provisional  time  levels  in  the  corrector  are  eliminated  by  substituting 
the  following  expansions: 

if1  -  F  (uf1)  .  F"  +  (f )“  (ufi  -  uj> 

+  J  <*4)  (uf1  -  U®)2  +  .  . . 
au  i  1  1 
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- -.1*-' _*_»  ,\j«  ‘»*  *_•  '.►  *„*  *.•  * _•■  ’.*  %»_\V  "-•.***  '.*  ',•  * *,> 


Fn+^  ■  F  (un+1)  ■  Fn  +  /JZ\n  /jtH+I  TTn  » 

i+1  *1+1/  Fi+1  +  Vi+l  (Ui+l  Ui+i) 


in  which 


+  I  (^T)"  (“JtJ  -  “?+i>2  +  ••• 

3U  i 


(B.4) 


or 

an  -  A  -  n  x  n  Jacobian  matrix 


-  3rd  order  n  x  n  x  n  tensor. 


However,  from  the  predictor. 


\vm 


9 


ui+1  ~  u°  *  -  0  <Fi  -  fJ.j) 


C  <Ci  -  r°> 


and  thus 


Fn+1  -  Fn  -  o(— )  (Fn  -  f”  )  +  2—  CFn  -  Fn  'i2  j. 

i  i  °  3U  A  lFi  Fi-1)  +  T  (^2)±  (Fi  Fi-i>  + 

—  n  0  7  n  (B.6) 

Fn+^  *  F11  -  a(— )  (Fn  -  Fn)  +  —  /JLZ.'i  /Fn  pn.2  . 

i+1  *i+i  q3Uj1+1  (Fi+l  V  +  T  CFi+l  “  Fi>  +  ••• 

Taylor  series  expansions  may  be  used  to  give, 

uf1  -  »?  ♦  at  .t»  +  +  T"  ”ttt’  +  ••• 

F"+l  -  F;  +  AX  Fx“  +  ^  Pxx"  +  ^  Fm»  +  ...  (,.7) 

’U  -  ¥  *°  +... 


(B.5) 


(B.6) 


(B.7) 


Ax  _  n 

-  p  4. 

6  xxx t  +  **• 


Substituting  Eqns.  (B.6),  (B.7)  into  (B.3)  will  give,  to  the  lowest 


order  remaining  terms, 
A  2 

r„  .  „  .Ax  ,  _ 


F  u  4  7  4.  Ax  r  p  i  At  Ax  0  e  9  /3Fv  v  * 

LUt  +  Fx+  6  {Fxxx}  3^(3^)  Fx} 

*  <  6  Uttt  +  Fx  (2  3U2  Fxx  +  I  Fx)} 

■  -  —  (u  -  Jl  fiZ  f  nin 
2lUtt  3x(3uV}Ji 


(B.8) 


However,  from  the  exact  equation,  (B.l) 


3U  9F 

at  "ax 


and  thus 
Utt 


—  (-F  ) 

at  '  x' 


a  ,3F  .  a  ,aF  . 

"  ax  (au  ut*  "  ax  (au  V 


(B.9) 


The  right  hand  side  of  Eqn.  (B.8)  is  then  identically  zero  and  the 
difference  algorithm  generates  an  approximation  of  the  form, 
ut  +  Fx  +  0(Ax2)  +  0(AxAt)  +  0(At2)  -  0 


(B.10) 


Appendix  C 


SPLITTING 

The  concept  o£  operator  factoring  or  tine  splitting  was  introduced 
in  Section  III  as  an  efficient  method  of  solving  multi-dimensional 
system  of  equations. 

It  remains  to  show  that  this  does  not  affect  the  formal  accuracy 
of  the  numerical  method.  Lomax ,  (Ref  58),  has  introduced  a  simple 
procedure  in  order  to  prove: 

(1)  A  symmetric  operator  sequence  is  a  necessary  condition  for 
second  order  temporal  accuracy, 

(2)  In  the  linear  case,  a  symmetric  operator  sequency  recovers 
the  full  second  order  accuracy  of  the  algorithm. 

In  what  follows,  the  spatial  accuracy  is  unaffected  by  the  time 
splitting  and  does  not  need  to  be  considered  in  detail.  The  two  dimen¬ 
sional  equation 


55L  +  il  +  lG  .  0 

3t  3x  +  3y  U 

may  be  written  in  quasi-linear  form  as 


(C.l) 


£  +  a  |H  +  a  |H.«  o 

3t  x  3x  y  3y 


(C.2) 


where 


A  3F 

Ax  “  3U 


A  -!£ 
y  3U 


If  the  Jacobian  matrices,  A  A  are  assumed  to  be  constant,  the  system 

x  y 


is  then  linear.  A  second  order  accurate  solution  to  Eqn.  (0.2)  must 


reproduce  the  terms  of  a  Taylor  series  expansion  up  to  order  two. 


U11*1  »  u"  +  At  +  *£_  ^  +  0(it3) 


(C.3) 


But  from  Eqn.  (C.2) 


Ut  "  <Ax  3x  +  Ay  V  U 


U  •  (A2  3  +  (A  A  +  A  A  )  3  +  A2  3  )U 

tt  x  xx  '  x  y  y  x  xy  y  yy 


(0.4) 


A  numerical  approximation  to  Eqn.  (C.3)  may  be  generated  by  time 


splitting  as  the  symmetric  sequence. 


U 


0+1  -  Ly( ix(At)  Ly( |^)  Un 


(C.5) 


or 


U 


v  * u 


L  u  ,  ir  - 

x 


L  U 
y 


** 


The  Ly  operator  solves  the  equation 


12.  +  A  12.  -  o 

3t  y  3y 


(C.6) 

by  the  MacCormack  predictor-corrector  to  give  the  second  order  accurate 


approximation 


U 


iy<f% 


where 


.  .At.  _  .  At  .  .  .  1  /At. 2  .2  . 

Ly  <2~>  -  1  +  r  A,  5,  +  2  <r>  \6yy 

The  lx  operator  similarly  solves 
3U  .  .  3U  n 

if  +  A*  s  ■ 0 


to  give 

** 

U 


-  L  (At)  U 
x 


where 


(At)  -  I  +  At  A„  6„  +  i  (At)2  A2  6. 


X  XX 

Then  Eqn.  (C.5)  may  be  written  as 


X  XX 


U 


n+1 


/At 


1  .At. 2  .2 


(C.7) 


(C.8) 


2.2, 


(I+(^)  A  6+4-  (^)‘  A  6  )(I  +  At  A  A  +  f  At^TA  ) 
2  y  y  2  2  y  yy  x  x  2  x  xx 


(C.9) 


+  (r>.  Ay  +  7  (j1)2  Ay  «yy)u" 

which  may  be  expanded  to  give 


U 


n+1 


(C.10) 


-  {  I  +  At  (A  A  +  A  6  ) 
y  y  xx 

+  *2-  <Ay  {yy  +  +  yjy,  +  +  0(it3) 
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Appendix  D 


STABILITY  CRITERIA  FOR  THE  TRANSFORMED 
NAVIER- STOKES  EQUATIONS 


The  purpose  of  this  appendix  is  to  show  how  the  stability  analysis 
of  MacCormack  and  Baldwin  (Ref  55)  may  be  generalized  for  the  trans¬ 
formed  equations  written  in  computational  variables.  The  dominant 
inviscid  limit  will  be  considered  in  some  detail  and  a  simple  general¬ 
ization  then  given  for  the  viscous  corrections.  Attention  will  be 
restricted  to  the  L ^  operator,  Eqn.  (3.36).  Inviscid  terms  in  Eqn. 
(3.37)  may  be  written  in  non-conservative  form  as: 

It  «>  +  <V  +{y  y  <yGx)  '°  (D-1) 

where 


p 

pu 

A 

pv 

pu 

F  - 

2 

pu 

G„  - 

puv 

pv 

I 

puv 

I 

2 

pv 

-PE. 

puE 

m  m 

pvE 

The  spectral  radius  of  the  amplification  matrix,  G,  associated  with 
Eqn.  (D.l) ,  may  be  more  conveniently  found  if  a  change  of  dependent 
variables  is  made,  (Ref  50:  Sec  13.4): 


V  - 


(D.2) 


Equation  (D.l)  can  then  be  written  in  quasi-llnear  form  as 


fj  <V)  +  A  H  +  Fj  V  -  0 


(D.3) 


where 


U  =  £  u  +  l  v 
c  x  y 


Following  Ricthmyer  (Ref  50?  See  8.4),  the  lower  order  term  F^V  may  be 
neglected.  The  amplification  matrix  for  Eqn.  (D.3)  may  be  determined 
as  outlined  in  Section  III  and  has  been  given  by  MacCormack  (Ref  54) 


as 

G  -  I  -  i  ||  A  sin  a  -  |  (||)2  {  A(l-e~ia)  }2  (D.4) 

where 

i  -  /-I  I  -  identity  matrix 

a  ■  k^A£  “  phase  angle 

To  find  the  stability  restriction,  one  must  find  the  largest  eigenvalue 
of  the  matrix  A,  which  now  differs  from  Reference  (55)  due  to  the  pres- 
sence  of  the  transformation  derivatives: 


-  U  ,  U  ,  U  +  /.  2.-2 
A  c’  c’  c  - 


(D.5) 


of  which 


(XA) 

A  max 


U  + 
c 


c 


(D.6) 


W 

'.•■V 


*  N  * 


The  eigenvalues  of  the  amplification  matrix  are  then 


l  +  i|fxA  sin  «,+!(§>  ni  (1-s-10)) 


-icu,2 


2  VAC' 


For  stability, 

IM  “  [1-(||-)2  (1-cos  a)]2  +  [(^|-)  *A  sin  “]2  1  1  (D. 7) 

Substituting  (^A)max  into  Eqn.  (D.7)  and  differentiating  with  respect 
to  a  to  find  the  extremal  values  gives 


|X(XA,a) 


a-ir/2 


|  X  |  <1 

'max  — 


(D.8) 


i* 


From  Eqn.  (D.8),  for  stability 


A L 


(Atr).  <  ^ - 

C  inv  —  (A.) 

A  max 


AS£ 


where 


AS, 


AC 


U 


A  2+c  2 

x  y 


.  uE 


/.  2  2 

»&x  +S 


(D.9) 


Equation  (D.9)  may  be  compared  with  results  given  by  MacCormack, 
(Ref  55). 

Ax 


(At  ) 


x  inv  u  +c 


(D.10) 


<:■ 

§i 

?v • 1 


*3 


The  similarity  in  form  suggests  that  the  viscous  corrections  may  be 


generalized  by  a  simple  change  of  variables: 

2 


(At  ) 


Ax 


becomes 


x'diff  2  yy/(Prp) 


AS' 


^Vdiff  ’  2  ,v  [  t  , 

*p?  +  p r)/p 


and 


<«*>,»d 


AxAy 

»-Ap/p 


(Ref  55) 


(Ref  55) 


(D.ll) 


(D.12) 


(D.13) 
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becomes 


<AVmd 


AS  AS 

s  n 


/| (y+5) (X+Xt) | h 


(D.14; 


where 


AS 


An 


/  2  2 
*  n  +n 
x  y 


Corresponding  expressions  for  the  L  operator  may  be  derived  in  a 


similar  manner. 


Appendix  E 

THROTTLE  RATIO  CORRECTION  TO  AREA  RATIO 
As  mentioned  In  Section  4.4,  the  throttle  ratio  settings  given 
in  Reference  1  are  based  upon  the  exit  area  of  the  holes  In  the  cowling. 
Slnde  this  is  not  the  physical  choking  area,  a  procedure  must  be  used 
to  determine  the  choking  area  if  the  numerical  simulation  is  to  be 
run  at  the  same  condition  as  the  experiment. 

Referring  to  Figure  35,  throttle  ratio,  TR,  is  based  upon  the 
exit  area 

A  ■  KirR.,h  (E.l) 

e  x 

where  k  is  a  constant  which  accounts  for  the  area  reduction  due  to  the 
presence  of  struts  at  the  outlet.  The  choking  area  is  calculated  as 
A*  -  irs  (Rx  +  R2)  (E.2) 


The  two  can  not  be  related  because  the  constant  k  was  not  reported  in 

Ref  1.  One-dimensional  relationships  may  be  used  to  determine  the 

throat  area,  such  that  the  pressure  values  agree  with  the  results 

of  Figure  7  from  the  experiment.  For  any  steady  state  flow,  the  mass 

fliix  across  all  crossectional  areas  of  the  duct  must  be  identically 

equal.  The  mass  flux,  captured  by  the  inlet  in  the  supercritical 

regime,  m  ,  was  calculated  by  applying  the  Navier-Stokes  equations 
csp 

only  to  the  upstream  region  and  may  be  considered  known. 

Since  the  area  ratio  and  pressure  at  station  7  were  given,  this 
station  may  be  used  to  infer  properties  at  the  downstream  throat.  A 


one-dimensional  area  ratio  may  be  written  as 


^  Ml(Y+1) 

_  cr  • 


1  r  2  ni  Y_1 

M^Ty+iT  (1  +  2  M7)J 


(E.3) 


where  A*  is  the  area  required  to  choke  the  throat  from  the  values 
eff 

specified  at  station  7.  If  there  is  little  total  pressure  loss  from 

station  7  to  the  throat,  p  ■  p  ,  A.  will  be  the  desired  area 

/  eff 

after  correcting  for  boundary  layer  displacement  effects.  The  mass 
flux  across  station  7  is  given  as 

*7- 'l  1=  m7  *  +  *r  m5  a7  ■  *c  (E-‘> 

»T 

c7 

with  Aj  *  irR^ 

P7  -  Pt  (P7/Pt  )  (from  Figure  7) 

Equation  (E.4)  may  be  solved  for  which  is,  in  turn,  substituted  into 

Eqn  (E.3)  to  give  A*  .  The  actual  geometric  area  will  be  somewhat 

eff 


due  to  the  presence  of  the  boundary  layers 


larger  than  Aa 

eff 


A*  ■ 

*  t 


eff 


(E.5) 


The  displacement  correction  was  obtained  from  preliminary  aft-diffuser 
calculations  in  which  mass-averaged  quantities  were  computed  from  the 
numerical  solution  and  substituted  into  Eqn  (E.4)  written  at  the  nozzle 
throat.  After  solving  for  A*  ,  with  the  geometric  area  known,  the 


eff 

displacement  correction  was  set  as  e,  ■  .969. 

a 

Now  that  the  two  areas  have  been  related,  the  blockage  constant, 

k,  can  be  computed  once  and  for  all.  For  TR  -  1.42 

A  -  TR  x  A  -  7.668  cm2 
e  c 

2 

where  A£  is  given  as  5.4  cm  .  To  match  the  resulting  pressure  at 

station  7  (Fig.  7)  requires 
2 

A^  -  6.367  cm  . 

Upon  using  Eqn  (E.2)  and  trignometric  relationships  from  Figure  (E.l) 
R^  ■  2  cm 
R2  ■  R^  -  s  cos  0 
s  ■  h  sin  0 

the  geometry  is  completely  determined.  Eqn  (E.l)  can  then  be  solved 
to  calculate  the  blockage  constant 
k  »  .655  B  .66 

Once  this  constant  has  been  fixed,  the  two  areas  are  directly  related 
by  geometry.  For  k  ■  .66,  TR  ■  1.42  corresponds  to  an  area  ratio 


AR  -  Ag/A*  -  1.162 

In  a  similar  manner,  TR  ■  ,97  corresponds  to  AR  •  .834. 


After  the  numerical  solution  has  been  obtained,  the  validity  of 
the  assumptions  made  in  the  analysis  may  be  checked  by  computing  one 
dimensional  mass-averaged  properties  from  the  full  solution.  For 


TR  -  1.42 


Pt  /Pt  *  -995 

*■*  i 


ed  "  A*  /A*  "  *967 
a  eff 
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SO.  ABSTRACT  (Continue  on  roreteo  elde  It  neceeeery  and  Identify  by  Block  number) 

-^The  unsteady,  compressible,  Reynolds-averaged  Navier-Stokes  equations 
were  solved  for  the  flow  field  about  an  external  compression  axi-symmetric 
inlet  with  a  length  to  diameter  ratio,  L/Dsl5.88,  at  Mach  2.0  and  a  Reynolds 
number  based  on  diameter,  Re_  ■  2.36x10®,  operating  in  the  near-critical  and 
subcritical  flow  regimes.  The  near-critical  solution  reached  a  stable  steady 
state  while  the  subcritical  solutions  attained  a*.  unstable  bounded  oscillatory 
state,  characterized  by  large  amplitude  pressure  oscillations  and  traveling 
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shock  waves.  This  phenomenon  is  a  result  of  a  shear  layer  instability  combined 
with  a  closed-loop  feedback  of  reflected  disturbances  and  the  naturally  occur¬ 
ring  self-sustained  oscillations  are  commonly  known  as  buzz.  Numerical  results 
are  given  in  terms  of  Mach  contours,  velocity  field  plots,  pressure-time  traces 
at  selected  stations,  as  well  as  mass  flux  and  other  mass-averaged  quantities 
along  the  duct  length.  Comparison  with  experiment  is  also  given. 
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